{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 第四周作业：对活动进行聚类\n",
    "\n",
    "第四周和第五周的作业均数据来源于Kaggle竞赛：Event Recommendation Engine Challenge，根据提供的信息判断用户对某个事件是否感兴趣。\n",
    "\n",
    "[竞赛官网](https://www.kaggle.com/c/event-recommendation-engine-challenge/data)\n",
    "\n",
    "### 样本描述\n",
    "活动描述信息在events.csv文件：共110维特征。\n",
    "\n",
    "前9列：event_id, user_id, start_time, city, state, zip, country, lat, and lng。\n",
    "\n",
    "- event_id：活动的id, \n",
    "- user_id：创建活动的用户的id .  \n",
    "- city, state, zip, and country： 活动地点 (如果知道的话).\n",
    "- lat and lng： floats（活动地点的经度和纬度）\n",
    "- start_time： 字符串，ISO-8601 UTC time，表示活动开始时间\n",
    "\n",
    "后101列为词频：count_1, count_2, ..., count_100，count_other。\n",
    "\n",
    "- count_N：活动描述出现第N个词的次数\n",
    "- count_other：除了最常用的100个词之外的其余词出现的次数\n",
    "\n",
    "### 作业要求：\n",
    "1. 根据活动的关键词（count_1, count_2, ..., count_100，count_other属性）做聚类，可采用KMeans聚类。\n",
    "2. 尝试K=10，20，30，..., 100, 并计算各自CH_scores。\n",
    "\n",
    "提示：由于样本数目较多，建议使用MiniBatchKMeans。\n",
    "\n",
    "文件说明：\n",
    "1.\t可以先运行0. EDA.ipynb，看一下竞赛所有数据的情况；\n",
    "2.\t总体活动的数目太多（300w+记录），可以只需对训练集train.csv和测试集test.cv出现的活动（13418条记录）举行聚类即可。运行1. Users_Events.ipynb可得到只在训练集train.csv和测试集test.cv出现的活动，可自己修改代码存为csv格式，在进行聚类。\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "#导入聚类和分析结果需要的包\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "\n",
    "from sklearn.decomposition import PCA\n",
    "from sklearn.cluster import MiniBatchKMeans\n",
    "from sklearn import metrics\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>user</th>\n",
       "      <th>event</th>\n",
       "      <th>invited</th>\n",
       "      <th>timestamp</th>\n",
       "      <th>interested</th>\n",
       "      <th>not_interested</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>3044012</td>\n",
       "      <td>1918771225</td>\n",
       "      <td>0</td>\n",
       "      <td>2012-10-02 15:53:05.754000+00:00</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>3044012</td>\n",
       "      <td>1502284248</td>\n",
       "      <td>0</td>\n",
       "      <td>2012-10-02 15:53:05.754000+00:00</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>3044012</td>\n",
       "      <td>2529072432</td>\n",
       "      <td>0</td>\n",
       "      <td>2012-10-02 15:53:05.754000+00:00</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>3044012</td>\n",
       "      <td>3072478280</td>\n",
       "      <td>0</td>\n",
       "      <td>2012-10-02 15:53:05.754000+00:00</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>3044012</td>\n",
       "      <td>1390707377</td>\n",
       "      <td>0</td>\n",
       "      <td>2012-10-02 15:53:05.754000+00:00</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "      user       event  invited                         timestamp  interested  \\\n",
       "0  3044012  1918771225        0  2012-10-02 15:53:05.754000+00:00           0   \n",
       "1  3044012  1502284248        0  2012-10-02 15:53:05.754000+00:00           0   \n",
       "2  3044012  2529072432        0  2012-10-02 15:53:05.754000+00:00           1   \n",
       "3  3044012  3072478280        0  2012-10-02 15:53:05.754000+00:00           0   \n",
       "4  3044012  1390707377        0  2012-10-02 15:53:05.754000+00:00           0   \n",
       "\n",
       "   not_interested  \n",
       "0               0  \n",
       "1               0  \n",
       "2               0  \n",
       "3               0  \n",
       "4               0  "
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#先读数据看看\n",
    "train=pd.read_csv('train.csv')\n",
    "train.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>user</th>\n",
       "      <th>event</th>\n",
       "      <th>invited</th>\n",
       "      <th>timestamp</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>1776192</td>\n",
       "      <td>2877501688</td>\n",
       "      <td>0</td>\n",
       "      <td>2012-11-30 11:39:01.230000+00:00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>1776192</td>\n",
       "      <td>3025444328</td>\n",
       "      <td>0</td>\n",
       "      <td>2012-11-30 11:39:01.230000+00:00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>1776192</td>\n",
       "      <td>4078218285</td>\n",
       "      <td>0</td>\n",
       "      <td>2012-11-30 11:39:01.230000+00:00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>1776192</td>\n",
       "      <td>1024025121</td>\n",
       "      <td>0</td>\n",
       "      <td>2012-11-30 11:39:01.230000+00:00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>1776192</td>\n",
       "      <td>2972428928</td>\n",
       "      <td>0</td>\n",
       "      <td>2012-11-30 11:39:21.985000+00:00</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "      user       event  invited                         timestamp\n",
       "0  1776192  2877501688        0  2012-11-30 11:39:01.230000+00:00\n",
       "1  1776192  3025444328        0  2012-11-30 11:39:01.230000+00:00\n",
       "2  1776192  4078218285        0  2012-11-30 11:39:01.230000+00:00\n",
       "3  1776192  1024025121        0  2012-11-30 11:39:01.230000+00:00\n",
       "4  1776192  2972428928        0  2012-11-30 11:39:21.985000+00:00"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "test=pd.read_csv('test.csv')\n",
    "test.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "发现train和test中只包含对应每个活动ID的用户ID。\n",
    "\n",
    "本次作业最终需要解决的问题是判断某用户是否会对某活动感兴趣，可能是因为对于整体的解决思路还不清晰，所以我觉得本次的作业要求比较模糊。\n",
    "\n",
    "由于这次作业只是对活动的关键词做聚类，再加上最终评价的CH指标是内部评价。所以我认为对于本次作业，只需获得train里的所有event，然后将event.csv中有相应的样本取出构成一个新样本集，最后对其101维活动关键词聚类即可。\n",
    "\n",
    "根据提示和文件说明运行了1. User_Events.ipynb，但报了不少错，所以先要对生成新样本集的代码进行修改。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 修改原User_Events.ipynb代码，从events.csv中获取train和test中的相应样本，并输出为新样本集events_train.csv"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "number of uniqueEvents :13418\n"
     ]
    }
   ],
   "source": [
    "#用set保存不重复的event\n",
    "uniqueEvents = set()\n",
    "\n",
    "for filename in [\"train.csv\", \"test.csv\"]:\n",
    "    f = open(filename, 'r')\n",
    "    \n",
    "    #忽略第一行（列名）\n",
    "    f.readline().strip().split(\",\")\n",
    "    \n",
    "    for line in f:    #对每条记录\n",
    "        cols = line.strip().split(\",\")\n",
    "        uniqueEvents.add(cols[1])   #第二列为活动ID\n",
    "        \n",
    "    f.close()\n",
    "\n",
    "n_uniqueEvents = len(uniqueEvents)\n",
    "\n",
    "print(\"number of uniqueEvents :%d\" % n_uniqueEvents)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "获得了如作业提示所述的事件数，现在遍历原始大数据集，用每条记录和set中的记录比较，将一致的信息全部保存到一个新的文件。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "13418\n",
      "totally cost 21.596235275268555\n"
     ]
    }
   ],
   "source": [
    "import time\n",
    "time_start=time.time()\n",
    "\n",
    "'''最简单的双重循环搜索原始数据集，速度慢到令人发指，等了1个小时还卡着，后来经过测试，这个方法大概需要十几个小时才能遍历完毕\n",
    "import linecache\n",
    "events = linecache.getlines(\"events.csv\")\n",
    "columns = events[0].strip().split(',')\n",
    "print(columns)\n",
    "del events[0]\n",
    "\n",
    "data1=[]\n",
    "\n",
    "for e in events:\n",
    "    e = e.strip().split(',')\n",
    "    for i in uniqueEvents:\n",
    "        if int(i) == int(e[0]):\n",
    "            data1.append(e)\n",
    "    \n",
    "target_data=pd.DataFrame(data=data1,columns=sample_columns)\n",
    "\n",
    "print(target_data.shape)\n",
    "\n",
    "print(type(i))\n",
    "print(type(e[0]))\n",
    "target_data.head()\n",
    "'''\n",
    "# 研究了很久后突然想起set有个in的基本操作，一下就把耗时降低到了几十秒\n",
    "events=list()\n",
    "# 读取原始数据集\n",
    "f = open(\"events.csv\", 'r')\n",
    "# 获得特征名称列表并跳过该行\n",
    "column = f.readline().strip().split(\",\")\n",
    "    \n",
    "for fLine in f:    #对每条记录\n",
    "    events_cols = tuple(fLine.strip().split(\",\"))\n",
    "    if events_cols[0] in uniqueEvents:#如果其event在train、test集中存在，就添加\n",
    "        events.append(events_cols)\n",
    "        \n",
    "'''使用with open耗时基本一样\n",
    "with open(\"events.csv\",\"r\") as f: \n",
    "    for fLine in f:\n",
    "        fLine = tuple(fLine.strip().split(','))\n",
    "        if fLine[0] in uniqueEvents:\n",
    "            events.add(fLine)\n",
    "'''\n",
    "print(len(events))\n",
    "time_end=time.time()\n",
    "print('totally cost',time_end-time_start)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "这一步研究了非常久的时间……按照作业的意思这个代码应该存在于1. User_Events.ipynb，但实际上不存在，而那个文件其实是研究用户和事件之间的关系，和本次作业没有直接联系。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_events=pd.DataFrame(data=events,columns=column)\n",
    "df_events.to_csv('train_events.csv',index=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "终于得到了所需要的train相关的events的dataframe，感动。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>event_id</th>\n",
       "      <th>user_id</th>\n",
       "      <th>start_time</th>\n",
       "      <th>city</th>\n",
       "      <th>state</th>\n",
       "      <th>zip</th>\n",
       "      <th>country</th>\n",
       "      <th>lat</th>\n",
       "      <th>lng</th>\n",
       "      <th>c_1</th>\n",
       "      <th>...</th>\n",
       "      <th>c_92</th>\n",
       "      <th>c_93</th>\n",
       "      <th>c_94</th>\n",
       "      <th>c_95</th>\n",
       "      <th>c_96</th>\n",
       "      <th>c_97</th>\n",
       "      <th>c_98</th>\n",
       "      <th>c_99</th>\n",
       "      <th>c_100</th>\n",
       "      <th>c_other</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>684921758</td>\n",
       "      <td>3647864012</td>\n",
       "      <td>2012-10-31T00:00:00.001Z</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>2</td>\n",
       "      <td>...</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>9</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>244999119</td>\n",
       "      <td>3476440521</td>\n",
       "      <td>2012-11-03T00:00:00.001Z</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>2</td>\n",
       "      <td>...</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>7</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>3928440935</td>\n",
       "      <td>517514445</td>\n",
       "      <td>2012-11-05T00:00:00.001Z</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>0</td>\n",
       "      <td>...</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>12</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>2582345152</td>\n",
       "      <td>781585781</td>\n",
       "      <td>2012-10-30T00:00:00.001Z</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>1</td>\n",
       "      <td>...</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>8</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>1051165850</td>\n",
       "      <td>1016098580</td>\n",
       "      <td>2012-09-27T00:00:00.001Z</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>1</td>\n",
       "      <td>...</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>9</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>5 rows × 110 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "     event_id     user_id                start_time city state  zip country  \\\n",
       "0   684921758  3647864012  2012-10-31T00:00:00.001Z  NaN   NaN  NaN     NaN   \n",
       "1   244999119  3476440521  2012-11-03T00:00:00.001Z  NaN   NaN  NaN     NaN   \n",
       "2  3928440935   517514445  2012-11-05T00:00:00.001Z  NaN   NaN  NaN     NaN   \n",
       "3  2582345152   781585781  2012-10-30T00:00:00.001Z  NaN   NaN  NaN     NaN   \n",
       "4  1051165850  1016098580  2012-09-27T00:00:00.001Z  NaN   NaN  NaN     NaN   \n",
       "\n",
       "   lat  lng  c_1   ...     c_92  c_93  c_94  c_95  c_96  c_97  c_98  c_99  \\\n",
       "0  NaN  NaN    2   ...        0     1     0     0     0     0     0     0   \n",
       "1  NaN  NaN    2   ...        0     0     0     0     0     0     0     0   \n",
       "2  NaN  NaN    0   ...        0     0     0     0     0     0     0     0   \n",
       "3  NaN  NaN    1   ...        0     0     0     0     0     0     0     0   \n",
       "4  NaN  NaN    1   ...        0     0     0     0     0     0     0     0   \n",
       "\n",
       "   c_100  c_other  \n",
       "0      0        9  \n",
       "1      0        7  \n",
       "2      0       12  \n",
       "3      0        8  \n",
       "4      0        9  \n",
       "\n",
       "[5 rows x 110 columns]"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df_events=pd.read_csv('train_events.csv')\n",
    "df_events.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>c_1</th>\n",
       "      <th>c_2</th>\n",
       "      <th>c_3</th>\n",
       "      <th>c_4</th>\n",
       "      <th>c_5</th>\n",
       "      <th>c_6</th>\n",
       "      <th>c_7</th>\n",
       "      <th>c_8</th>\n",
       "      <th>c_9</th>\n",
       "      <th>c_10</th>\n",
       "      <th>...</th>\n",
       "      <th>c_92</th>\n",
       "      <th>c_93</th>\n",
       "      <th>c_94</th>\n",
       "      <th>c_95</th>\n",
       "      <th>c_96</th>\n",
       "      <th>c_97</th>\n",
       "      <th>c_98</th>\n",
       "      <th>c_99</th>\n",
       "      <th>c_100</th>\n",
       "      <th>c_other</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>2</td>\n",
       "      <td>0</td>\n",
       "      <td>2</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>...</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>9</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>2</td>\n",
       "      <td>0</td>\n",
       "      <td>2</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>...</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>7</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>...</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>12</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>2</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>...</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>8</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>2</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>...</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>9</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>5 rows × 101 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "   c_1  c_2  c_3  c_4  c_5  c_6  c_7  c_8  c_9  c_10   ...     c_92  c_93  \\\n",
       "0    2    0    2    0    0    0    0    0    0     0   ...        0     1   \n",
       "1    2    0    2    0    0    0    0    0    0     0   ...        0     0   \n",
       "2    0    0    0    0    0    0    0    0    0     0   ...        0     0   \n",
       "3    1    0    2    1    0    0    0    0    0     0   ...        0     0   \n",
       "4    1    1    0    0    0    0    0    2    0     0   ...        0     0   \n",
       "\n",
       "   c_94  c_95  c_96  c_97  c_98  c_99  c_100  c_other  \n",
       "0     0     0     0     0     0     0      0        9  \n",
       "1     0     0     0     0     0     0      0        7  \n",
       "2     0     0     0     0     0     0      0       12  \n",
       "3     0     0     0     0     0     0      0        8  \n",
       "4     0     0     0     0     0     0      0        9  \n",
       "\n",
       "[5 rows x 101 columns]"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#先删除无用特征\n",
    "train=df_events.drop(columns=['event_id','user_id','start_time','city','state','zip','country','lat','lng'])\n",
    "train.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "有101个特征要进行聚类"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 一个参数点（聚类数据为K）的模型\n",
    "def K_cluster_analysis(K, train):\n",
    "    start = time.time()\n",
    "    \n",
    "    #K-means,在训练集上训练\n",
    "    mb_kmeans = MiniBatchKMeans(n_clusters = K)\n",
    "    mb_kmeans.fit(train)\n",
    "    \n",
    "    # 作业要求使用Calinski-Harabasz Index评价，不过我也想看一下Silhouette Coefficient\n",
    "    chi = metrics.calinski_harabaz_score(train,mb_kmeans.predict(train))\n",
    "    sc = metrics.silhouette_score(train,mb_kmeans.predict(train))\n",
    "    \n",
    "    \n",
    "    end = time.time()\n",
    "    print(\"K:{},CH_score: {}, Silhouette_Coefficient: {},time elaps:{}\".format(K, chi, sc, int(end-start)))\n",
    "    \n",
    "    return chi,sc"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "K:10,CH_score: 1341.2486266720525, Silhouette_Coefficient: 0.4500938937481291,time elaps:11\n",
      "K:20,CH_score: 527.7504389151641, Silhouette_Coefficient: 0.2468259414674305,time elaps:9\n",
      "K:30,CH_score: 361.02803979379496, Silhouette_Coefficient: 0.21947212239669361,time elaps:9\n",
      "K:40,CH_score: 686.9296704836412, Silhouette_Coefficient: 0.15894549409451572,time elaps:8\n",
      "K:50,CH_score: 261.5039743888655, Silhouette_Coefficient: 0.1436901382161007,time elaps:8\n",
      "K:60,CH_score: 200.47355261571084, Silhouette_Coefficient: 0.10659341040006358,time elaps:8\n",
      "K:70,CH_score: 192.08474513260526, Silhouette_Coefficient: 0.110989713667812,time elaps:8\n",
      "K:80,CH_score: 221.57369389245954, Silhouette_Coefficient: 0.10599506405019879,time elaps:8\n",
      "K:90,CH_score: 162.48382324728897, Silhouette_Coefficient: 0.08712390702254362,time elaps:8\n",
      "K:100,CH_score: 115.11871737224466, Silhouette_Coefficient: 0.0810774328370229,time elaps:9\n"
     ]
    }
   ],
   "source": [
    "# 设置超参数（聚类数目K）范围\n",
    "Ks = range(10,101,10)\n",
    "\n",
    "chis = []\n",
    "scs = []\n",
    "for K in Ks:\n",
    "    chi, sc = K_cluster_analysis(K, train)\n",
    "    chis.append(chi)\n",
    "    scs.append(sc)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x5226668>]"
      ]
     },
     "execution_count": 27,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYAAAAD8CAYAAAB+UHOxAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzt3XmYVNWd//H3F5odka1FZBGiiMYtYodFcKFUBFwzMYlOEok/IjNx1zEKro8ao46ZqETFoBh1YlzGMZFExBBAUbbQiAugBtygRaUNiwoKIt/549z+dQHddNO13Oq6n9fz1FNVp05XfamnqE/dc+4919wdERFJniZxFyAiIvFQAIiIJJQCQEQkoRQAIiIJpQAQEUkoBYCISELVGQBm9oCZrTazxTU8dpmZuZl1ju6bmY03s+Vm9pqZ9UvrO8rMlkWXUdn9Z4iIyK6qzxbAg8Dw7RvNrAdwPLAirXkE0Ce6jAEmRH07AtcBA4D+wHVm1iGTwkVEJDN1BoC7zwLW1PDQ7cDlQPqRZKcCD3swD2hvZl2BE4Bp7r7G3dcC06ghVEREJH9KGvJHZnYK8IG7v2pm6Q91A1am3a+I2mpr36nOnTt7r169GlKiiEhiLVy48BN3L62r3y4HgJm1Bq4ChtX0cA1tvpP2mp5/DGH4iJ49e1JeXr6rJYqIJJqZvV+ffg3ZC2gfoDfwqpm9B3QHXjazPQm/7Huk9e0OrNpJ+w7cfaK7l7l7WWlpnQEmIiINtMsB4O6vu/se7t7L3XsRvtz7uftHwGTgrGhvoIHAenf/EHgOGGZmHaLJ32FRm4iIxKQ+u4E+CswF+ppZhZmN3kn3KcA7wHLgPuBcAHdfA9wILIguN0RtIiISEyvk5aDLyspccwAiIrvGzBa6e1ld/XQksIhIQikAREQSSgEgIpJQRRkAa9bADTfAokVxVyIiUrgadCRwoWvaFK6/HrZsgcMOi7saEZHCVJRbALvvDmVlMGNG3JWIiBSuogwAgFQK5s+Hzz+PuxIRkcJU1AGwZQu89FLclYiIFKaiDYDBg6FZMw0DiYjUpmgDoHVrGDRIASAiUpuiDQAIw0Avvwxr18ZdiYhI4Sn6AHCHF16IuxIRkcJT1AEwYAC0aqVhIBGRmhR1ADRvDkceqQAQEalJUQcAhGGgJUvg44/jrkREpLAkIgAAZs6Mtw4RkUJT9AFw2GFhaQgNA4mIbKvoA6CkBI4+WgEgIrK9og8ACMNAb78N778fdyUiIoUjMQEAmgcQEUmXiAA48EAoLdUwkIhIukQEQJMmMHRoCAD3uKsRESkMiQgACMNAH3wAy5bFXYmISGGoMwDM7AEzW21mi9PabjOzN83sNTP7o5m1T3tsnJktN7O3zOyEtPbhUdtyMxub/X/KzlXNA2gYSEQkqM8WwIPA8O3apgEHufshwD+AcQBm9k3gDODA6G/uMbOmZtYUuBsYAXwTODPqmzf77gvduysARESq1BkA7j4LWLNd21/dfUt0dx7QPbp9KvCYu29y93eB5UD/6LLc3d9x983AY1HfvDELWwEzZ8LWrfl8ZRGRwpSNOYD/Bzwb3e4GrEx7rCJqq609r1Ip+OQTWLy47r4iIsUuowAws6uALcAjVU01dPOdtNf0nGPMrNzMyisrKzMpbwdDh4ZrDQOJiGQQAGY2CjgJ+KH7/9+5sgLokdatO7BqJ+07cPeJ7l7m7mWlpaUNLa9GPXuGuQAFgIhIAwPAzIYDVwCnuPvGtIcmA2eYWQsz6w30Af4OLAD6mFlvM2tOmCienFnpDZNKhTOEbdlSd18RkWJWn91AHwXmAn3NrMLMRgN3AbsB08zsFTO7F8DdlwBPAEuBqcB57v51NGF8PvAc8AbwRNQ371Ip+PTTcK5gEZEkK6mrg7ufWUPzpJ30vwm4qYb2KcCUXaouB9LnAfr3j7cWEZE4JeZI4Cp77AEHH6x5ABGRxAUAhGGgl16CTZvirkREJD6JDYAvvoD58+OuREQkPokMgKOOCiuEahhIRJIskQHQvj0cfrgCQESSLZEBAGEYaN482LAh7kpEROKR6AD46iuYPTvuSkRE4pHYABg8GJo10zCQiCRXYgOgTRsYOFABICLJldgAgDAMtHAhrFsXdyUiIvmX+ADYuhVmzYq7EhGR/Et0AAwYAK1aaRhIRJIp0QHQogUMGaIAEJFkSnQAQBgGev11WL067kpERPJLAZAK188/H2sZIiJ5l/gA6NcP2rXTMJCIJE/iA6CkBI4+WgEgIsmT+ACAMAy0bBmsXBl3JSIi+aMAoHoeYObMeOsQEcknBQBw0EHQubOGgUQkWRQAhJPDDB0aAsA97mpERPJDARBJpcIcwNtvx12JiEh+KAAiVfMAGgYSkaRQAET69IFu3RQAIpIcdQaAmT1gZqvNbHFaW0czm2Zmy6LrDlG7mdl4M1tuZq+ZWb+0vxkV9V9mZqNy889pOLOwFaB5ABFJivpsATwIDN+ubSww3d37ANOj+wAjgD7RZQwwAUJgANcBA4D+wHVVoVFIUimorIQlS+KuREQk9+oMAHefBazZrvlU4KHo9kPAaWntD3swD2hvZl2BE4Bp7r7G3dcC09gxVGI3dGi41jCQiCRBQ+cAurj7hwDR9R5Rezcg/XjaiqittvYdmNkYMys3s/LKysoGltcwe+8N++yjABCRZMj2JLDV0OY7ad+x0X2iu5e5e1lpaWlWi6uPVCqsDPr113l/aRGRvGpoAHwcDe0QXVetpl8B9Ejr1x1YtZP2gpNKwfr1sGhR3JWIiORWQwNgMlC1J88o4Om09rOivYEGAuujIaLngGFm1iGa/B0WtRUczQOISFLUZzfQR4G5QF8zqzCz0cAtwPFmtgw4ProPMAV4B1gO3AecC+Dua4AbgQXR5YaoreB06QIHHqgAEJHiV1JXB3c/s5aHjq2hrwPn1fI8DwAP7FJ1MUmlYNIk2LwZmjePuxoRkdzQkcA1SKVg40b4+9/jrkREJHcUADU4+uhwZLCGgUSkmCkAatChQzhXsAJARIqZAqAWqRTMnRuGgkREipECoBapVJgEnjMn7kpERHJDAVCLIUOgpETDQCJSvBQAtWjbFgYMUACISPFSAOxEKgULFoSlIUREio0CYCdSKdi6FV58Me5KRESyTwGwEwMHQsuWGgYSkeKkANiJli1h8GAFgIgUJwVAHVIpePVV+OSTuCsREckuBUAdUqlw/fzzsZYhIpJ1CoA6lJXBbrtpGEhEio8CoA4lJXDUUQoAESk+CoB6SKXgrbfggw/irkREJHsUAPVQNQ8wc2a8dYiIZJMCoB4OOQQ6dtQwkIgUFwVAPTRpEk4WP306uMddjYhIdigA6imVghUr4J134q5ERCQ7FAD1VDUPoGEgESkWCoB66tsXunZVAIhI8VAA1JNZ2AqYMUPzACJSHDIKADO7xMyWmNliM3vUzFqaWW8zm29my8zscTNrHvVtEd1fHj3eKxv/gHxKpWD1ali6NO5KREQy1+AAMLNuwIVAmbsfBDQFzgBuBW539z7AWmB09CejgbXuvi9we9SvUdE8gIgUk0yHgEqAVmZWArQGPgRSwJPR4w8Bp0W3T43uEz1+rJlZhq+fV716Qe/eCgARKQ4NDgB3/wD4FbCC8MW/HlgIrHP3LVG3CqBbdLsbsDL62y1R/04Nff24pFJhZdCvv467EhGRzGQyBNSB8Ku+N7AX0AYYUUPXqinTmn7t7zCdamZjzKzczMorKysbWl7OpFKwbh288krclYiIZCaTIaDjgHfdvdLdvwKeAo4A2kdDQgDdgVXR7QqgB0D0+O7Amu2f1N0nunuZu5eVlpZmUF5uDB0arjUMJCKNXSYBsAIYaGato7H8Y4GlwEzg9KjPKODp6Pbk6D7R4zPcG98OlV27wgEHKABEpPHLZA5gPmEy92Xg9ei5JgJXAJea2XLCGP+k6E8mAZ2i9kuBsRnUHatUCl58ETZvjrsSEZGGs0L+EV5WVubl5eVxl7GDp56C734XXnopnDReRKSQmNlCdy+rq5+OBG6AY44JRwZrGEhEGjMFQAN07AiHHaYAEJHGTQHQQKkUzJkDX3wRdyUiIg2jAGigVCpMAs+ZE3clIiINowBooCFDoKREw0Ai0ngpABpot92gf38FgIg0XgqADKRSsGABfPpp3JWIiOw6BUAGUqmwKNyLL8ZdiYjIrlMAZGDQIGjRQsNAItI4KQAy0LJlOBJYASAijZECIEOpVFga+p//jLsSEZFdowDIUNVpIp9/PtYyGo0CXnpKJHEUABkqK4O2bTUMVB+XXBLWUdLZ1EQKgwIgQ82awVFHKQDq8te/wh13wKxZ8OSTdfcXkdxTAGRBKgVvvgmrVtXdN4k2bIB/+zfo2xf23x9uvBG2bo27KhFRAGRB1TzAzJnx1lGorrkG3nsP7rsv3F6yBP74x7irEhEFQBYceih06KBhoJr8/e9w553ws5/BkUfCD34A++0HN9ygrQCRuCkAsqBJk3CyeAXAtjZvhtGjw3mUb7kltDVtCldfDa+9BpMnx1ufSNIpALIklQrDHO++G3clhePWW2HxYpgwAdq1q24/80zYZ5+wFaDdQkXiowDIkqp5AG0FBG+8Ab/4RRjyOfnkbR8rKYGrroJFi+CZZ+KpT0QUAFmz//6w554KAAhj++ecE46PGD++5j4/+hH07q2tAJE4KQCyxCxsBcyYoS+0CRNg9my4/XbYY4+a+zRrBldeGZbTfu65/NYnIoECIItSKfjoo3BMQFKtWAFjx8KwYfDjH++871lnQc+ecP31Ck2ROGQUAGbW3syeNLM3zewNMxtkZh3NbJqZLYuuO0R9zczGm9lyM3vNzPpl559QOJI+D+AedvfcuhV++9uwVbQzzZvDuHEwbx787W/5qVFEqmW6BXAnMNXd9wcOBd4AxgLT3b0PMD26DzAC6BNdxgATMnztgtO7N/TqldwAeOwxmDIFbropvA/1cfbZ0L27tgJE4tDgADCzdsBRwCQAd9/s7uuAU4GHom4PAadFt08FHvZgHtDezLo2uPIClUqFI4KTdpDTJ5/AhRfCgAFwwQX1/7sWLcKQ0ezZWlFVJN8y2QL4BlAJ/M7MFpnZ/WbWBuji7h8CRNdV04DdgJVpf18RtRWVVArWroVXX427kvy65BJYtw7uvz8c7LUrqg4Wu/763NQmIjXLJABKgH7ABHc/DNhA9XBPTWoaEd5ho9/MxphZuZmVV1ZWZlBePIYODddJGgaaOhV+//swnn/QQbv+9y1bwhVXwAsvhIuI5EcmAVABVLj7/Oj+k4RA+LhqaCe6Xp3Wv0fa33cHdlg/090nunuZu5eVlpZmUF489torHBOQlAD4/POw0ucBB4SDuxrqnHOgS5ewUqiI5EeDA8DdPwJWmlnfqOlYYCkwGRgVtY0Cno5uTwbOivYGGgisrxoqKjapVFj3/quv4q4k9666ClauDEM/LVo0/Hlat4af/xymTw/zASKSe5nuBXQB8IiZvQZ8C/glcAtwvJktA46P7gNMAd4BlgP3Aedm+NoFK5UKv4zLy+OuJLfmzoXf/AbOOw+OOCLz5/v3f4fS0nB0sIjkXkkmf+zurwBlNTx0bA19HTgvk9drLI45JlzPmAGDBsVaSs5s2gQ//WnYhfOXv8zOc7ZpA5ddFuYD5s2DgQOz87wiUjMdCZwDnTrBt75V3PMAN98MS5fCvffCbrtl73nPPTe8f5oLEMk9BUCOpFJhLPvLL+OuJPuWLAm/+v/1X2HkyOw+d9u2cOml4YCyYh9CE4mbAiBHUqkwTDJ3btyVZNfXX4ehn3btwknec+H888MZ1rQVIJJbCoAcOfLIcEBUsQ0D3X13GJ+/884wYZsL7dqFA8smTw7nDBCR3FAA5Ei7dvDtbxdXALz/fljCecSIMPyTSxdcALvvrq0AkVxSAORQKhVOiv7ZZ3FXkjn3cMAXhPX+61rpM1Pt28NFF8Ef/xjOHywi2acAyKGRI2HLFjjlFFi9uu7+heyRR8KJW26+GfbeOz+vefHFYQ+jX/wiP68nkjQKgBwaPBgefjiMmZeVNd69Wiorw5fxoEFhN8186dAhrDD65JNhzyMRyS4FQI79+Mdhd1AzGDIEHnww7op23UUXwaefNmylz0xdcklYJkJbASLZpwDIg379YOHCEABnnx2WTti8Oe6q6ueZZ+DRR8OaP9/8Zv5fv1OnsFvo448n+1SbIrmgAMiTzp3DssmXXQb33BMmiD8s8KXwPvssrM9z4IFhqee4/Md/QKtW4UxjIpI9CoA8KimB224Lv6gXLYLDDy/sA8XGjYMPPghDP82bx1dHaWmYe/jDH+Af/4ivDpFiowCIwRlnhC/+Vq3g6KPDCdQL7Xy4s2eHLZULLyyMRdkuuyyEULYWnhMRBUBsDjkEFiyAY48NwyznnFM46wZ9+WVY7qFnz8KZfO3SJbxPv/89vP123NWIFAcFQIw6doS//CVMsE6aFLYGKiririqMtb/5Ztgyads27mqq/fznYRjt5pvjrkSkOCgAYta0afiV/dRTYXnlww8PZxOLy+uvwy23hN1XTzghvjpqstdeMGYMPPQQvPde3NWINH4KgALxne+EZSM6dAjDQuPH539eoGqlzw4d4Pbb8/va9XX55dCkibYCRLJBAVBADjgA5s8PS0hcdBGMGgVffJG/1x8/PoTQ+PFh//tC1L07jB4Nv/sdrFgRdzUijZsCoMDsvntYAO2GG8KE5+DB+RnuePdduPpqOPFE+MEPcv96mRg7Nlzfemu8dYg0dgqAAtSkCVxzDfz5z/DOO2EdoenTc/d6VSt9Nm2an5U+M9WzZzii+v77w3EKItIwCoACduKJYVfRLl1g2DD41a9yMy/w8MMwbVqY/O3RI/vPnwvjxsHWrdoKEMmEAqDA9ekT5gX+5V/CbpBnngkbNmTv+T/+OCy4NmRI2M++sejVC846CyZOLPwlNUQKlQKgEWjbFp54IvxC/5//CUfmLl+enee+8MIQKPfdF4aeGpMrrwznW7jttrgrEWmcMv4vb2ZNzWyRmf0lut/bzOab2TIze9zMmkftLaL7y6PHe2X62kliBldcERaUW7UqnG7y2Wcze87Jk0OwXHMN7L9/durMp332gR/9CO69N2zJiMiuycZvvouAN9Lu3wrc7u59gLXA6Kh9NLDW3fcFbo/6yS46/vhwYpm99w5zBDfdFMbCd9X69fCzn8HBB4d96xurK6+ETZvgv/4r7kpEGp+MAsDMugMnAvdH9w1IAU9GXR4CTotunxrdJ3r82Ki/7KLevWHOnHBi9quvhu9+N5ywZVeMHQsffRSWoIhzpc9M7bdfmBe5++5w5jIRqb9MtwDuAC4Hqn6DdgLWufuW6H4F0C263Q1YCRA9vj7qLw3QujX893+HI3b//GcYMADeeqt+fztrVhg2ufjiMJTU2F11VThg7te/jrsSkcalwQFgZicBq919YXpzDV29Ho+lP+8YMys3s/JK/aTbKbPwJf63v8E//xm+zJ9+eud/8+WXYeXR3r3DwWbF4IAD4Pvfh7vuCu+DiNRPJlsAg4FTzOw94DHC0M8dQHszK4n6dAdWRbcrgB4A0eO7A2u2f1J3n+juZe5eVlpamkF5yXHMMeGUk337wmmnwbXX1j4vcOON4aQqEydCmzZ5LTOnrrkGPv8c7rgj7kpEGo8GB4C7j3P37u7eCzgDmOHuPwRmAqdH3UYBVb9JJ0f3iR6f4V5op0FpvHr0gBdfDEfI3ngjnHIKrFu3bZ9XX4X//E/4yU/guONiKTNnDjwQTj89rGO0dm3c1Yg0DrnY8/sK4FIzW04Y458UtU8COkXtlwJjc/DaidayZZjUveceeO65MCS0ZEl4bMuWsIhap07Fu8fM1VeHyfA774y7EpHGwQr5R3hZWZmXl5fHXUajNHt2+EX82Wfw4INhQbmf/zzs9/+978VdXe585zswcya8/35YWE8kicxsobuX1dWvkR37KfU1eHCYFzjkkPCFP25cGBY6/fS6/7Yxu/bacIzDb34TdyUihU8BUMT22guefx7OPTfcvueewl/pM1OHHQYnnxx2Cd3VYyNEkkYBUOSaNw8HSb33HnTrVmf3onDNNWEi+O67465EpLApABKi2H/5p/v2t2HEiDDZ/fnncVcjUrgUAFKUrr02HBQ2YULclYgULgWAFKWBA8NJdG67DTZujLsakcKkAJCide21YYG4e++NuxKRwqQAkKI1eDCkUuHo5y++iLsakcKjAJCidt114WQx990XdyUihUcBIEXtqKPg6KPDyeO//DLuakQKiwJAit6114bTaE6aVHdfkSRRAEjRGzo0zAfccks4faSIBAoAKXpmYSugoiIsjCcigQJAEuH448OxAb/8JWzeHHc1IoVBASCJULUVsGJFOJeyiCgAJEGGD4eyMrjpJvjqq7irEYmfAkASo2or4N13w8FhWihOkk4BIIly0kkwZEg4fWSnTnDCCeEUksuXx12ZSP4pACRRzGD69HA5//wwJ3DxxdCnD+y3H1xyCUybpt1FJRl0TmBJvHfegSlT4JlnwvmEN22CNm3guOPgxBPDuQW6d4+7SpH6q+85gRUAImk2boQZM6oDYcWK0H7ooTByZAiEAQOgpCTeOkV2RgEgkiF3WLo0BMGUKfDSS/D119ChQ9ijaOTIcN25c9yVimxLASCSZevWhfmBKVPCZfXqMKcwYEDYMhg5MpyUPkmn35TCpAAQyaGtW+Hll6u3DhYsCFsMXbuGOYMTTwxzCO3axV1pYdm4Ed5+O+x1tf1lzZqwh9bw4eGy334K04bKeQCYWQ/gYWBPYCsw0d3vNLOOwONAL+A94PvuvtbMDLgTGAlsBH7i7i/v7DUUANJYrF4NU6eGQHjuOVi/Hpo1gyOPrJ476Ns3GV9on35a+5f8qlXb9u3cGfbdN1zatg2T8G+9FR7r1as6DFIp2G23vP9TGq18BEBXoKu7v2xmuwELgdOAnwBr3P0WMxsLdHD3K8xsJHABIQAGAHe6+4CdvYYCQBqjLVtgzpzqieTFi0N7797VQ0XHHAOtWsVaZkbWrq35C3758hCG6fbcs/pLPv2yzz7Qvv2Oz/3uuyFEp04Nu+t+/nkI0/Stg4MPTkaYNlTeh4DM7GngruhyjLt/GIXE8+7e18x+G91+NOr/VlW/2p5TASDFYMWK6nmD6dPDMEirVuFX7T77QOvW1Zc2bWq/n367VSto2jR3NbvDJ5/U/iW/Zs22/Xv02PaLPf1227YNr2PzZpg9O4TB1Knw2muhvWvX6jA47jjo2LHhr1GM8hoAZtYLmAUcBKxw9/Zpj6119w5m9hfgFnd/KWqfDlzh7uXbPdcYYAxAz549D3///fczrk+kUHz5JbzwQgiD556Djz4KgdCQtYlattw2LOoKkNruN2sGK1fu+CX/6afVr9WkCey9d82/5Hv3zt/WzKpV1VsHf/1rmJhv0iRMxI8YEQLh8MNDW5LlLQDMrC3wAnCTuz9lZutqCYBngJu3C4DL3X1hbc+tLQBJiq++Cieu37gRNmwI11WXuu7X92/qCpmSkvBlXtOXfK9e0Lx5Xt6KetuyJUy+P/tsCITy8rDl0rkzDBsWwmDYMOjSJe5K86++AZDR4Sxm1gz4X+ARd38qav7YzLqmDQFVjQhWAD3S/rw7sN2UkEgyNWsWLrncayg9ZNJDYtOmcKRzz56N6wC3khIYNChcbrgBKivDbrpTp4athD/8IfTr1y+EwYgR4ZwQjenfmGuZTAIb8BBhwvfitPbbgH+mTQJ3dPfLzexE4HyqJ4HHu3v/nb2GtgBEpCG2boVXXqmeO5gzJxzEt/vuYc5g+PCwEGCPHnU/V2OUj72AhgAvAq8TdgMFuBKYDzwB9ARWAN9z9zVRYNwFDCfsBnr29uP/21MAiEg2rFsXJuCrAqGiIrQfeGD1ZPKRR0KLFvHWmS06EExEpAZVS3xUhcGsWWFvo9atw55Zxx8PgwfDIYeEYbnGSAEgIlIPGzbA88+HyeRnnw2rw0IIhP794YgjqucaOnWKtdR6UwCIiDTAihUwd26YN5g7FxYtCnscQTia+4gjqkPhgAMKc5dTBYCISBZs3Bh2Ma0KhDlzwkFyEI5kHjiwOhT69y+MJSsUACIiOeAeDpSbM6f6smRJaG/SJCxTURUIRxwRjq3I97IVCgARkTxZvx7mz68OhHnz4LPPwmN77LFtIBx+eDiKO5fyciCYiIiE4wuGDQsXCMccLF267VbCn/4UHmvWLITAoEHVobDXXvHUrS0AEZE8WL06bBlUBcKCBWFtKAjrLKUHQqa7oGoISESkgG3eHI5WrppYnj0bPvggPNa6NZx8Mjz2WMOeW0NAIiIFrHnzsNdQ//5w0UWhbeXK6kBo0yb3NSgAREQKRI8e4fL97+fn9QrwEAYREckHBYCISEIpAEREEkoBICKSUAoAEZGEUgCIiCSUAkBEJKEUACIiCVXQS0GYWSXwftx1ZKgz8EncRRQQvR/b0vtRTe/FtjJ5P/Z299K6OhV0ABQDMyuvz5ocSaH3Y1t6P6rpvdhWPt4PDQGJiCSUAkBEJKEUALk3Me4CCozej23p/aim92JbOX8/NAcgIpJQ2gIQEUkoBUAWmVkPM5tpZm+Y2RIzuyhq72hm08xsWXTdIe5a88XMmprZIjP7S3S/t5nNj96Lx82sedw15ouZtTezJ83szegzMijhn41Lov8ni83sUTNrmaTPh5k9YGarzWxxWluNnwcLxpvZcjN7zcz6ZaMGBUB2bQH+w90PAAYC55nZN4GxwHR37wNMj+4nxUXAG2n3bwVuj96LtcDoWKqKx53AVHffHziU8L4k8rNhZt2AC4Eydz8IaAqcQbI+Hw8Cw7drq+3zMALoE13GABOyUYACIIvc/UN3fzm6/RnhP3g34FTgoajbQ8Bp8VSYX2bWHTgRuD+6b0AKeDLqkqT3oh1wFDAJwN03u/s6EvrZiJQArcysBGgNfEiCPh/uPgtYs11zbZ+HU4GHPZgHtDezrpnWoADIETPrBRwGzAe6uPuHEEIC2CO+yvLqDuByYGt0vxOwzt23RPcrCAGZBN8AKoHfRUNi95vg1qemAAAB2klEQVRZGxL62XD3D4BfASsIX/zrgYUk9/NRpbbPQzdgZVq/rLw3CoAcMLO2wP8CF7v7p3HXEwczOwlY7e4L05tr6JqU3dBKgH7ABHc/DNhAQoZ7ahKNbZ8K9Ab2AtoQhjm2l5TPR11y8n9HAZBlZtaM8OX/iLs/FTV/XLW5Fl2vjqu+PBoMnGJm7wGPETbt7yBsupZEfboDq+IpL+8qgAp3nx/df5IQCEn8bAAcB7zr7pXu/hXwFHAEyf18VKnt81AB9Ejrl5X3RgGQRdEY9yTgDXf/ddpDk4FR0e1RwNP5ri3f3H2cu3d3916Eyb0Z7v5DYCZwetQtEe8FgLt/BKw0s75R07HAUhL42YisAAaaWevo/03V+5HIz0ea2j4Pk4Gzor2BBgLrq4aKMqEDwbLIzIYALwKvUz3ufSVhHuAJoCfhg/89d99+8qdomdkxwGXufpKZfYOwRdARWAT8yN03xVlfvpjZtwgT4s2Bd4CzCT/CEvnZMLPrgR8Q9p5bBPyUMK6diM+HmT0KHENY9fNj4DrgT9TweYhC8i7CXkMbgbPdvTzjGhQAIiLJpCEgEZGEUgCIiCSUAkBEJKEUACIiCaUAEBFJKAWAiEhCKQBERBJKASAiklD/B3XMbdcZlRUCAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "#用图看一下\n",
    "plt.plot(Ks, np.array(chis), 'b-')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x528da20>]"
      ]
     },
     "execution_count": 28,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX0AAAD8CAYAAACb4nSYAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzt3Xl8VPW9//HXJwkhJOwSFsOSICEUxYJEsFoDCgqyJb319oLllrb2qq0Lre3jd7X6uPZna2+v9z5cWrWVasX+enGpWlYVKxWDtQJBcAMDAVnCImERBAJk+fz+mAEHCGQgQ04y834+HvNIzjnfc+Yzw+F9Tr5z5nvM3RERkcSQFHQBIiLSeBT6IiIJRKEvIpJAFPoiIglEoS8ikkAU+iIiCUShLyKSQBT6IiIJRKEvIpJAUoIu4HidOnXy7OzsoMsQEWlWli1btsPdM+tr1+RCPzs7m5KSkqDLEBFpVsxsQzTt1L0jIpJAFPoiIglEoS8ikkAU+iIiCUShLyKSQKIKfTMbbWalZlZmZnecot21ZuZmlh+ezjazSjNbEX78LlaFi4jI6av3kk0zSwYeBa4CyoGlZjbb3Vce164NcBuw+LhNrHX3gTGqV0REGiCaM/0hQJm7r3P3w8CzQGEd7X4O3A8cjGF9UdtVuYt737yX97a9F8TTi4g0C9GEfhawKWK6PDzvKDMbBPRw97l1rJ9jZsvN7E0zu7yuJzCzG8ysxMxKKioqoq39GEmWxM+Lf84zHz5zRuuLiCSCaELf6ph39G7qZpYEPAj8uI52W4Ge7j4IuB2YYWZtT9iY+zR3z3f3/MzMer9FXKf2ae0Znj2cWaWzzmh9EZFEEE3olwM9Iqa7A1siptsAFwALzWw9cAkw28zy3f2Qu+8EcPdlwFqgbywKr0tRXhEf7/iY0h2lZ+spRESatWhCfymQa2Y5ZpYKTARmH1no7nvcvZO7Z7t7NvAOMMHdS8wsM/xBMGbWG8gF1sX8VYRNyJsAoLN9EZGTqDf03b0auAWYD6wCnnf3j8zsXjObUM/qBcD7ZvYe8AJwk7vvamjRJ9OjXQ8GdxvMzI9nnq2nEBFp1qIaZdPdXwZePm7ef5yk7fCI318EXmxAfaetMK+Qexbew7Z92+jaumtjPrWISJMXd9/ILepXhOPMKZ0TdCkiIk1O3IX+BZ0vIKd9jvr1RUTqEHehb2YU9Svi9XWvs+/wvqDLERFpUuIu9CHUr3+o5hDzy+YHXYqISJMSl6F/Wc/LOKfVOcws1VU8IiKR4jL0U5JSGNd3HPNWz6OqpirockREmoy4DH0IXcWz++BuFm1cFHQpIiJNRtyG/lW9ryItJY1ZH+sqHhGRI+I29DNSM7j6vKuZWToTd69/BRGRBBC3oQ+hq3g27tnIim0rgi5FRKRJiOvQH993PEmWpC9qiYiExXXoZ2ZkcmmPSzUAm4hIWFyHPoTG2H/v0/dY/9n6oEsREQlc3Id+Yb/Q7Xx1FY+ISAKEfp+OfTg/83z164uIkAChD6EvahVvKGZX5Vm7f4uISLOQEKFfmFdIjdcwb/W8oEsREQlUVKFvZqPNrNTMyszsjlO0u9bM3MzyI+bdGV6v1MxGxaLo0zX43MFktcnSAGwikvDqDf3wjc0fBa4B+gOTzKx/He3aALcBiyPm9Sd0I/XzgdHAY0dulN6YkiyJCXkTmF82n8qqysZ+ehGRJiOaM/0hQJm7r3P3w8CzQGEd7X4O3A8cjJhXCDzr7ofc/ROgLLy9RlfUr4j9VftZ8MmCIJ5eRKRJiCb0s4BNEdPl4XlHmdkgoIe7zz3ddRvL8OzhtG3ZVpduikhCiyb0rY55R0cwM7Mk4EHgx6e7bsQ2bjCzEjMrqaioiKKk05eanMqY3DHMXj2bmtqas/IcIiJNXTShXw70iJjuDmyJmG4DXAAsNLP1wCXA7PCHufWtC4C7T3P3fHfPz8zMPL1XcBoK8wrZvn87izcvrr+xiEgciib0lwK5ZpZjZqmEPpidfWShu+9x907unu3u2cA7wAR3Lwm3m2hmLc0sB8gFlsT8VUTpmj7X0CKphcbiEZGEVW/ou3s1cAswH1gFPO/uH5nZvWY2oZ51PwKeB1YCrwI3u3tgfSvt0tpxRc4VzPxYY+yLSGKyphZ++fn5XlJScta2/9ulv+UHL/+AlT9YyZcyv3TWnkdEpDGZ2TJ3z6+vXUJ8IzfShLzQHycai0dEElHChX5W2ywuPvdi9euLSEJKuNCH0FU8izcvZuvnW4MuRUSkUSVk6Bf1KwJgdunselqKiMSXhAz9/pn9Oa/DeerXF5GEk5Chb2YU9StiwScL+PzQ50GXIyLSaBIy9CHUr3+45jCvlr0adCkiIo0mYUP/0h6X0im9k8bYF5GEkrChn5yUzPi+45m3eh5VNVVBlyMi0igSNvQhdBXPnkN7eHPDm0GXIiLSKBI69Ef2HkmrlFYaY19EEkZCh356i3RG9RnFrNJZGoBNRBJCQoc+hK7i2bR3E8u3LQ+6FBGRsy7hQ39c33EkWZLG4hGRhJDwod8pvRNf7flVhb6IJISED32AorwiPtj+Aet2rwu6FBGRs0qhDxT2KwTQVTwiEvcU+kDvDr0Z0HmABmATkbgXVeib2WgzKzWzMjO7o47lN5nZB2a2wszeMrP+4fnZZlYZnr/CzH4X6xcQK4V5hSzauIgdB3YEXYqIyFlTb+ibWTLwKHAN0B+YdCTUI8xw9wHuPhC4H3ggYtladx8YftwUq8JjrahfEbVey7zV84IuRUTkrInmTH8IUObu69z9MPAsUBjZwN33RkxmAM3um04XdbuI7m27awA2EYlr0YR+FrApYro8PO8YZnazma0ldKZ/W8SiHDNbbmZvmtnlDar2LDIzCvMKmV82nwNVB4IuR0TkrIgm9K2OeSecybv7o+5+HvDvwN3h2VuBnu4+CLgdmGFmbU94ArMbzKzEzEoqKiqirz7GCvMKqayu5PV1rwdWg4jI2RRN6JcDPSKmuwNbTtH+WaAIwN0PufvO8O/LgLVA3+NXcPdp7p7v7vmZmZnR1h5zw7KH0a5lO126KSJxK5rQXwrkmlmOmaUCE4Fj7ihuZrkRk2OBNeH5meEPgjGz3kAu0GS/AZWanMqY3DHMWT2HmtqaoMsREYm5ekPf3auBW4D5wCrgeXf/yMzuNbMJ4Wa3mNlHZraCUDfOlPD8AuB9M3sPeAG4yd13xfxVxFBRvyIqDlTwj/J/BF2KiEjMpUTTyN1fBl4+bt5/RPw+9STrvQi82JACG9voPqNpkdSCmR/P5Ks9vxp0OSIiMaVv5B6nbcu2jOg9gpkfz9QY+yISdxT6dSjMK2Tt7rWsrFgZdCkiIjGl0K/DhLzQRxUai0dE4o1Cvw7ntjmXIVlDNMa+iMQdhf5JFOUVsXTLUjbv3Rx0KSIiMaPQP4kjY+zPLp1dT0sRkeZDoX8SX+r0JXI75qpfX0TiikL/JI4MwPa3T/7GnoN7gi5HRCQmFPqnUNSviKraKl4tezXoUkREYkKhfwqXdL+EzPRMjbEvInFDoX8KyUnJTMibwMtrXuZwzeGgyxERaTCFfj0K8wrZe2gvC9cvDLoUEZEGU+jXY2TvkaS3SNcY+yISFxT69WjVohWjzhvFrNJZGoBNRJo9hX4UivoVsfnzzSzbuizoUkREGkShH4WxuWNJtmSNxSMizZ5CPwrnpJ/D5b0u17dzRaTZU+hHqTCvkA+3f0jZrrKgSxEROWNRhb6ZjTazUjMrM7M76lh+k5l9YGYrzOwtM+sfsezO8HqlZjYqlsU3psK80ABsuopHRJqzekPfzJKBR4FrgP7ApMhQD5vh7gPcfSBwP/BAeN3+wETgfGA08Fh4e81OToccLuxyobp4RKRZi+ZMfwhQ5u7r3P0w8CxQGNnA3fdGTGYAR65tLASedfdD7v4JUBbeXrNUlFfE3zf9nYr9FUGXIiJyRqIJ/SxgU8R0eXjeMczsZjNbS+hM/7bTWbe5KOxXSK3XMnf13KBLERE5I9GEvtUx74RvKbn7o+5+HvDvwN2ns66Z3WBmJWZWUlHRdM+iB3UdRI+2PTQAm4g0W9GEfjnQI2K6O7DlFO2fBYpOZ113n+bu+e6en5mZGUVJwTgyxv5f1/6VA1UHgi5HROS0RRP6S4FcM8sxs1RCH8wecw9BM8uNmBwLrAn/PhuYaGYtzSwHyAWWNLzs4BT1K6KyupLX1r4WdCkiIqet3tB392rgFmA+sAp43t0/MrN7zWxCuNktZvaRma0AbgemhNf9CHgeWAm8Ctzs7jVn4XU0moJeBbRPa6+reESkWbKmNohYfn6+l5SUBF3GKU1+aTKvlr3Ktp9sIyUpJehyREQws2Xunl9fO30j9wwU5hWys3Inb296O+hSREROi0L/DIzuM5rU5FQNwCYizY5C/wy0admGETkjNMa+iDQ7Cv0zVNSviHW71/Hh9g+DLkVEJGoK/TM0vu94AF3FIyLNikL/DHVr041Lul+ifn0RaVYU+g1QmFfIsq3L2LRnU/2NRUSaAIV+AxT1C402Mbt0dj0tRUSaBoV+A/Tr1I++5/RVv76INBsK/QYqyivijfVv8NnBz4IuRUSkXgr9BirsV0h1bTWvrHkl6FJEROql0G+goVlD6ZLRRWPsi0izoNBvoOSkZMb3Hc8ra17hUPWhoMsRETklhX4MFPUr4vPDn/PG+jeCLkVE5JQU+jEwovcIMlpkMOtjXcUjIk2bQj8G0lLSGN1nNLNKZ1HrtUGXIyJyUgr9GCnqV8TWfVsp2dK0bwAjIolNoR8jY3LHkGzJGotHRJq0qELfzEabWamZlZnZHXUsv93MVprZ+2a2wMx6RSyrMbMV4UfcjlfQsVVHhmUP07dzRaRJqzf0zSwZeBS4BugPTDKz/sc1Ww7ku/uFwAvA/RHLKt19YPgxgThWmFfIyoqVrNm5JuhSRETqFM2Z/hCgzN3Xufth4FmgMLKBu7/h7gfCk+8A3WNbZvNQmBd6W3S2LyJNVTShnwVEjh1cHp53MtcDkWMSpJlZiZm9Y2ZFda1gZjeE25RUVFREUVLT1Kt9LwZ2Hah+fRFpsqIJfatjXp03hjWzyUA+8N8Rs3u6ez5wHfCQmZ13wsbcp7l7vrvnZ2ZmRlFS01WUV8Tbm97m9XWv6/65ItLkRBP65UCPiOnuwJbjG5nZSOAuYIK7Hx2PwN23hH+uAxYCgxpQb5P3r1/+VzIzMrnq/13FwMcH8uS7T1JZVRl0WSIiQHShvxTINbMcM0sFJgLHXIVjZoOAxwkF/vaI+R3MrGX4907AZcDKWBXfFPXu0Jv1U9fzxPgncHe+N+d79HiwB3ctuIvNezcHXZ6IJDiLpgvCzMYADwHJwB/c/T4zuxcocffZZvY6MADYGl5lo7tPMLNLCR0MagkdYB5y9ydP9Vz5+fleUhIfX3BydxauX8jDix9mdulskpOSubb/tUwdOpVLul8SdHkiEkfMbFm4K/3U7Zpav3M8hX6kdbvX8ciSR3hy+ZPsPbSXIVlDmDp0Ktf2v5bU5NSgyxORZi7a0Nc3chtJ7w69eWDUA2y+fTOPXPMInx38jG++9E2yH8rmF8W/YPv+7fVvRESkgXSmH5Bar2V+2XweXvww89fOp2VySyYNmMTUoVMZ2HVg0OWJSDOj7p1mZFXFKn6z5Dc8/d7THKg6QEGvAqYOnUphXiHJSclBlycizYBCvxnaXbmbJ5c/ySNLHmHDng30ateLW4bcwvWDrqdDqw5BlyciTZhCvxmrrq1mdulsHl78MMUbiklvkc6UL0/htqG30a9Tv6DLE5EmSKEfJ1ZsW8GvF/+aGR/M4FDNIUadN4qpQ6cyqs8okkyfw4tIiEI/zmzfv51py6bx2NLH2LpvK33P6cutQ27l2wO/TevU1kGXJyIB0yWbcaZzRmfuLrib9T9cz//+0//SPq09t75yK1kPZHH7/Nv5ZPcnQZcoIs2AzvSbsXfK3+HhxQ/zwsoXqKmtYULeBKYOncrw7OGY1TVOnojEK3XvJJDNezfz2NLHeHzZ4+ys3MmFXS7ktiG3cd2A62jVolXQ5YlII1D3TgLJapvFfSPuY9OPNh0z0Nt5vz6P97a9F3R5ItKEKPTjSKsWrbj+out576b3WPCtBSQnJXPF01ewZPOSoEsTkSZCoR+HzIwrc65k0XcW0T6tPSP/OJJFGxYFXZaINAEK/TiW3T6bRd9ZxLltzmXUn0bx+rrXgy5JRAKm0I9zWW2zePPbb9KnYx/GzRjH3NVzgy5JRAKk0E8AXVp3YeG3FzKgywC+9tzX+PNHfw66JBEJiEI/QXRs1ZHX//V1hmYNZeKLE/nje38MuiQRCYBCP4G0S2vH/MnzuSL7CqbMnMLjJY8HXZKINLKoQt/MRptZqZmVmdkddSy/3cxWmtn7ZrbAzHpFLJtiZmvCjymxLF5OX0ZqBnOvm8vY3LHcNO8mHvzHg0GXJCKNqN7QN7Nk4FHgGqA/MMnM+h/XbDmQ7+4XAi8A94fX7QjcAwwFhgD3mJkGhg9YWkoaL/3LS1zb/1puf+127iu+L+iSRKSRpETRZghQ5u7rAMzsWaAQWHmkgbu/EdH+HWBy+PdRwF/dfVd43b8Co4FnGl66NERqcirPfP0ZWqW04u437mZ/1X7uu/I+jdkjEueiCf0sYFPEdDmhM/eTuR545RTrZh2/gpndANwA0LNnzyhKklhISUphetF0WqW04j/f+k/2H97PQ6MfUvCLxLFoQr+uBKhzlDYzmwzkA8NOZ113nwZMg9CAa1HUJDGSZEn8btzvSG+RzkOLH6KyupLfjv2t7s0rEqeiCf1yoEfEdHdgy/GNzGwkcBcwzN0PRaw7/Lh1F55JoXL2mBkPjHqAjNQM7lt0HweqDjC9aDopSdHsHiLSnETzv3opkGtmOcBmYCJwXWQDMxsEPA6MdvftEYvmA7+M+PD2auDOBlctMWdm/OLKX5DeIp27/nYXB6sPMuPrM0hNTg26NBGJoXqv3nH3auAWQgG+Cnje3T8ys3vNbEK42X8DrYE/m9kKM5sdXncX8HNCB46lwL1HPtSVpumnl/+UB0c9yIurXuRrz32NyqrKoEsSkRjSTVSkTr9f9ntunHsjV+RcwayJs3QfXpEmTjdRkQb5t8H/xh+/9kfeXP8mo/40ij0H9wRdkojEgEJfTmryhZN57trnWLp5KSP+OIKdB3YGXZKINJBCX07p6/2/zsyJM/lw+4cMf3o4n+77NOiSRKQBFPpSrzG5Y5h33TzW7V5HwfQCyveWB12SiJwhhb5EZUTvEbw2+TW27dtGwVMFfLL7k6BLEpEzoNCXqF3W8zIWfGsBew7t4fKnLqd0R2nQJYnIaVLoy2nJPzefhVMWUlVbRcH0At7/9P2gSxKR06DQl9M2oMsAir9dTIukFgyfPpySLfpehUhzodCXM5LXKY9F31lEu7R2jPjjCP6+8e9BlyQiUVDoyxnL6ZDDou8somvrrlz9p6tZsG5B0CWJSD0U+tIg3dt2p/jbxfTu0JuxM8Yyb/W8oEsSkVNQ6EuDdWndhYVTFnJB5wv42nNf48WVLwZdkoichEJfYuKc9HNY8K0FXJx1Md944Rv86f0/BV2SiNRBoS8x0y6tHfMnz2d49nC+9ZdvMW3ZtKBLEpHjKPQlplqntmbupLlck3sNN869kYffeTjokkQkgkJfYq5Vi1b85V/+wj996Z/44fwfMvWVqby18S0OVR+qf2UROat0ExU5a6prq7lxzo38YcUfAGiZ3JKh3YdS0LOAgl4FfKXHV3RzFpEYifYmKlGFvpmNBh4GkoEn3P1Xxy0vAB4CLgQmuvsLEctqgA/CkxvdfQKnoNCPPzsP7OStjW9RvKGY4o3FvLv1XWq9lmRLZvC5g48eBC7reRkdW3UMulyRZilmoW9mycBq4CqgnNC9bie5+8qINtlAW+AnwOzjQn+fu0d9OqfQj3+fH/qcf5T/I3QQ2FDM4s2LOVxzGIABnQdQ0Ct0ELi85+V0a9Mt4GpFmodoQz8lim0NAcrcfV14w88ChcDR0Hf39eFltWdUrSSUNi3bcPV5V3P1eVcDcLD6IEs2L2HRhkUUbyxm+orpPLr0UQByO+YePQAU9Cogu302ZhZk+SLNWjShnwVsipguB4aexnOkmVkJUA38yt1nnsa6kgDSUtKOnt3fxV1U11azfOvyo91BL616iSeXPwmEvgFc0KvgaJdQv079dBAQOQ3RhH5d/6NO59Pfnu6+xcx6A38zsw/cfe0xT2B2A3ADQM+ePU9j0xKPUpJSuDjrYi7OupgfX/pjar2Wj7Z/xKKNiyjeUMwbn7zBjA9mANApvdMxfwl8ucuXSU5KDvgViDRd0YR+OdAjYro7sCXaJ3D3LeGf68xsITAIWHtcm2nANAj16Ue7bUkMSZbEgC4DGNBlAD+4+Ae4O2t3r6V4Q/HRA8FLq14CoG3LtlzW47Kjfznkn5tPanJqwK9ApOmIJvSXArlmlgNsBiYC10WzcTPrABxw90Nm1gm4DLj/TIsVATAz+nTsQ5+OffjuoO8CUL63PPSZQLhL6M4FdwKhrqNLul9CQc8CvnH+Nzi/8/lBli4SuGgv2RxD6JLMZOAP7n6fmd0LlLj7bDO7GPgL0AE4CGxz9/PN7FLgcaCW0BfBHnL3J0/1XLp6R2Jhx4EdX1wmuqGY5duWA/D9/O9z7xX36tJQiTsxvU6/MSn05WzYeWAnP1v4Mx4reYwOaR2478r7+N5F31P/v8SNaENfwzBIQjgn/Rx+M+Y3LL9xORd0voCb5t3Exb+/WHf8koSj0JeEcmGXC3ljyhs8d+1zVByo4KtPfZXJL01m897NQZcm0igU+pJwzIxvnP8NPr75Y+6+/G5eWPkCeY/k8V9v/ZcGhZO4p9CXhJWRmsHPr/w5K29eycjeI7ljwR0M+O0AXl7zctCliZw1Cn1JeL079GbmxJm8+s1XSbIkxs4Yy7gZ4yjbVRZ0aSIxp9AXCRvVZxTvf/99/ueq/6F4QzHnP3Y+d75+J/sO7wu6NJGYUeiLREhNTuXHl/6Y0ltKmXTBJH7191+R90geMz6YQVO7vFnkTCj0RerQrU03phdN5+3vvk231t345kvfpGB6ASu2rQi6NJEGUeiLnMJXenyFJf+2hCfGP0HpjlIGTxvM9+d+n50HdgZdmsgZUeiL1CPJkrj+outZfetqbh1yK79/9/fk/iaXx5Y+Rk1tTdDliZwWhb5IlNqnteeh0Q+x4qYVDOw6kJtfvpnB0wZTvKE46NJEoqbQFzlNF3S+gAXfWsCf//nP7D64m2HThzHpxUmU7y0PujSRein0Rc6AmXFt/2tZdfMq7hl2DzM/nkneI3n8ctEvOVh9MOjyRE5KoS/SAOkt0vnZ8J+x6uZVjO4zmrv+dhcXPHYBc0rn6BJPaZIU+iIxkN0+mxe/8SKvTX6N1ORUJjw7gbEzxrJ65+qgSxM5hsbTF4mxqpoqHl36KPcsvIfKqkp+dMmPuLvgbtq0bHPWnvNwzWF2V+5mV+Wuo4/dB4+djpy/5+Ae2rZsS6f0TnRK70RmeubR3zuldyIz44vpDmkddN+BZkA3UREJ2Kf7PuXOBXfy1Iqn6Nq6K/ePvJ/JF07GzOps7+7sr9p/QnhHE+KnGirCMDq06kDHVh3pkBb62bZlWz4//DkV+yvYcWAHOw7sYH/V/jrXT7IkOrbqeOxB4fiDRMR0ZkYmGS0yTvo65exQ6Is0EUs2L+HWV25lyeYlXNrjUoacO+SkAV5VW3XS7aQmp9KxVcdjHkdC/FTz26W1I8nq78mtrKo8egA48qg4UHHK6era6jq31TK55Yl/NbQ68a+IrDZZ9D2nrw4QMRDT0Dez0cDDhO6R+4S7/+q45QWE7qF7ITDR3V+IWDYFuDs8+Qt3f/pUz6XQl3hU67U8veJp7n7jbj4/9HkomFtFBHZaxxPnHRfi6S3Sm1Q4ujt7Du059qCw/9QHid0Hd5+wnR5tezCu7zjG9x3PFTlXkJaSFsCraf5iFvpmlgysBq4CyoGlwCR3XxnRJhtoC/wEmH0k9M2sI1AC5AMOLAMGu/uJ//JhCn2R+FVVU8Wuyl1HDwJrdq3h5TUv89ra19hftZ/0Fulc1fsqxvUdx9jcsXRr0y3okpuNaEM/JYptDQHK3H1deMPPAoXA0dB39/XhZbXHrTsK+Ku77wov/yswGngmiucVkTjTIrkFXVp3oUvrLgAMyx7G9y76HoeqD7Fw/ULmrJ7DnNVzmFU6C4D8c/MZ33c84/qOY1DXQU3qL53mKppLNrOATRHT5eF50WjIuiKSIFqmtGRUn1E8MuYR1k9dz/s3vc8vr/wlqcmp/Gzhzxg8bTDdH+zOjXNuZO7quRyoOhB0yc1WNGf6dR1ao/30N6p1zewG4AaAnj17RrlpEYlHZsaALgMY0GUAd15+JxX7K3h5zcvMXTOXZz58hmnvTiMtJY2RvUcyLncc4/qOI6utziWjFU3olwM9Iqa7A1ui3H45MPy4dRce38jdpwHTINSnH+W2RSQBZGZkMmXgFKYMnMLhmsMUbyhmTmmoG2ju6rkwDwZ1HXS0G2jwuYOjulopUUXzQW4KoQ9yRwCbCX2Qe527f1RH2+nA3OM+yF0GXBRu8i6hD3J3nez59EGuiETD3Vm1YxVzV89lzuo5vL3pbWq9lq6tuzI2dyzj+45nZO+RZKRmBF1qo4j1JZtjCF2SmQz8wd3vM7N7gRJ3n21mFwN/AToAB4Ft7n5+eN3vAj8Nb+o+d3/qVM+l0BeRM7HzwE5eKXuFuavn8mrZq+w5tIeWyS25MudKxvUNdQP1bBe/3cf6cpaIJKyqmire2vjW0auBynaVAXBhlwuPdgMNyRoSV91ACn0RkbDSHaVHPwN4a+Nb1HgNnTM6MyZ3TOhLYdlX0D6tfbO+JFShLyJSh12Vu5hfNp85q+fwStkrfHbwMyA0THa31t3o2ror3dp0o2tG+GfrrnRr3e3o75npmU1yADqFvohIPaprq/n7xr+zePNitu3bxtZ9W0M/Pw/93HNozwnrJFkSnTM6f3GAiDxQHDed3iK90V5LLL+RKyISl1KSUhiWPYwb6OxdAAAEh0lEQVRh2cPqXH6g6gCf7vv0hIPB0el9W3nv0/f4dN+n1HjNCeu3SW1z4sEg4q+GI/POST+n0T5fUOiLiJxEeot0cjrkkNMh55Ttampr2Fm588SDwudb2bY/9PPdre+ydd/WOofBTklKoUtGFy7vdTnPfP3sjlKj0BcRaaDkpGQ6Z3Smc0ZnvsyXT9l23+F9bNu3rc6/HLq27nrWa1Xoi4g0otaprenTsQ99OvYJ5Pnj5yJVERGpl0JfRCSBKPRFRBKIQl9EJIEo9EVEEohCX0QkgSj0RUQSiEJfRCSBNLkB18ysAtgQdB0N1AnYEXQRTYjej2Pp/fiC3otjNeT96OXumfU1anKhHw/MrCSa0e4Shd6PY+n9+ILei2M1xvuh7h0RkQSi0BcRSSAK/bNjWtAFNDF6P46l9+MLei+OddbfD/Xpi4gkEJ3pi4gkEIV+A5lZDzN7w8xWmdlHZjY1PL+jmf3VzNaEf3YIutbGYmbJZrbczOaGp3PMbHH4vXjOzFKDrrGxmFl7M3vBzD4O7yNfSfB940fh/ycfmtkzZpaWSPuHmf3BzLab2YcR8+rcHyzk12ZWZmbvm9lFsahBod9w1cCP3f1LwCXAzWbWH7gDWODuucCC8HSimAqsipj+L+DB8HuxG7g+kKqC8TDwqrv3A75M6H1JyH3DzLKA24B8d78ASAYmklj7x3Rg9HHzTrY/XAPkhh83AL+NSQXurkcMH8As4CqgFOgWntcNKA26tkZ6/d3DO+6VwFzACH3ZJCW8/CvA/KDrbKT3oi3wCeHPziLmJ+q+kQVsAjoSumvfXGBUou0fQDbwYX37A/A4MKmudg156Ew/hswsGxgELAa6uPtWgPDPzsFV1qgeAv4PUBuePgf4zN2rw9PlhP7zJ4LeQAXwVLi76wkzyyBB9w133wz8D7AR2ArsAZaRuPvHESfbH44cJI+IyXuj0I8RM2sNvAj80N33Bl1PEMxsHLDd3ZdFzq6jaaJcMpYCXAT81t0HAftJkK6cuoT7qguBHOBcIINQF8bxEmX/qM9Z+b+j0I8BM2tBKPD/191fCs/+1My6hZd3A7YHVV8jugyYYGbrgWcJdfE8BLQ3s5Rwm+7AlmDKa3TlQLm7Lw5Pv0DoIJCI+wbASOATd69w9yrgJeBSEnf/OOJk+0M50COiXUzeG4V+A5mZAU8Cq9z9gYhFs4Ep4d+nEOrrj2vufqe7d3f3bEIf0P3N3b8JvAFcG26WEO8FgLtvAzaZWV541ghgJQm4b4RtBC4xs/Tw/5sj70dC7h8RTrY/zAa+Fb6K5xJgz5FuoIbQl7MayMy+CiwCPuCLfuyfEurXfx7oSWhn/2d33xVIkQEws+HAT9x9nJn1JnTm3xFYDkx290NB1tdYzGwg8ASQCqwDvkPoZCsh9w0z+7/AvxC66m058D1C/dQJsX+Y2TPAcEKjaX4K3APMpI79IXxgfITQ1T4HgO+4e0mDa1Doi4gkDnXviIgkEIW+iEgCUeiLiCQQhb6ISAJR6IuIJBCFvohIAlHoi4gkEIW+iEgC+f/WCzXg/2fvnwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(Ks, np.array(scs), 'g-')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "虽然CHI在60的时候有个小高峰，不过显而易见，在10左右是最佳的，所以在10这边进一步微调。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "K:2,CH_score: 3622.8916101544387, Silhouette_Coefficient: 0.7562667331685791,time elaps:11\n",
      "K:4,CH_score: 2470.153137006816, Silhouette_Coefficient: 0.5625773507337796,time elaps:10\n",
      "K:6,CH_score: 2093.0098236346325, Silhouette_Coefficient: 0.5048541746153251,time elaps:10\n",
      "K:8,CH_score: 1649.3859216574667, Silhouette_Coefficient: 0.47201625064329833,time elaps:10\n",
      "K:10,CH_score: 2945.882405612799, Silhouette_Coefficient: 0.4413896735798426,time elaps:10\n",
      "K:12,CH_score: 976.9162806341591, Silhouette_Coefficient: 0.3904577792475824,time elaps:10\n",
      "K:14,CH_score: 799.4724562295695, Silhouette_Coefficient: 0.35727557509813007,time elaps:10\n",
      "K:16,CH_score: 582.9881073001947, Silhouette_Coefficient: 0.23764899001053896,time elaps:9\n",
      "K:18,CH_score: 493.6511787256354, Silhouette_Coefficient: 0.26827105217779107,time elaps:9\n"
     ]
    }
   ],
   "source": [
    "# 设置超参数（聚类数目K）范围\n",
    "Ks2 = range(2,20,2)\n",
    "\n",
    "chis = []\n",
    "scs = []\n",
    "for K in Ks2:\n",
    "    chi, sc = K_cluster_analysis(K, train)\n",
    "    chis.append(chi)\n",
    "    scs.append(sc)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x50b8eb8>]"
      ]
     },
     "execution_count": 30,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYAAAAD8CAYAAAB+UHOxAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzt3XecVNX9//HXh6qCCuiKhBKIosEYRV0RLFHBgoKzmmL0F5VEDYmx9xZjiSaiYosGJdFIjIpdUbAQsCRR0MWCIihEjK6grIK9Ap/fH5+7X1ZY2DYzd2bn/Xw85jGzZ+/MfFaXee8959xzzN0REZHS0yrtAkREJB0KABGREqUAEBEpUQoAEZESpQAQESlRCgARkRKlABARKVEKABGREqUAEBEpUW3SLmBNNtxwQ+/du3faZYiIFJUZM2a85+5l9R1X0AHQu3dvKisr0y5DRKSomNn/GnKcuoBEREqUAkBEpEQpAERESpQCQESkRCkARERKlAJARKREKQBEREpUiwyADz+Es8+GuXPTrkREpHAV9IVgTfX553DllfDGG3DLLWlXIyJSmFrkGcDGG8Pxx8Ntt8HMmWlXIyJSmFpkAACceiqstx6cc07alYiIFKYWGwCdO0cITJgA06alXY2ISOGpNwDMbC0ze8bMXjSzWWZ2ftJ+k5nNN7MXklv/pN3M7Gozm2dmM81s21qvNcLM5ia3Ebn7scLxx0NZGfz2t7l+JxGR4tOQM4AvgcHuvjXQHxhqZgOT753q7v2T2wtJ2z5A3+Q2EhgDYGZdgHOBHYABwLlm1jl7P8qqOnaM2UBTpsRNRERWqDcAPHySfNk2ufkanlIB/D153jSgk5l1A/YGJrv7YndfAkwGhjav/Pr96lfQs2cEga+pahGREtOgMQAza21mLwCLiA/x6cm3Lkq6ea4ws/ZJW3fgrVpPr0raVteeU2utBb/7HUyfDg88kOt3ExEpHg0KAHdf5u79gR7AADPbEjgT+C6wPdAFOD053Op6iTW0f4OZjTSzSjOrrK6ubkh59RoxAvr2jbGA5cuz8pIiIkWvUbOA3P0D4HFgqLsvTLp5vgT+RvTrQ/xl37PW03oAC9bQvvJ7jHX3cncvLyurd0ezBmnbFi64AF56CW6/PSsvKSJS9BoyC6jMzDolj9cG9gDmJP36mJkB+wMvJ0+ZAByWzAYaCHzo7guBR4C9zKxzMvi7V9KWFwceCFttFd1BX3+dr3cVESlcDTkD6AY8ZmYzgWeJMYAHgVvM7CXgJWBD4MLk+EnA68A84C/AbwDcfTHw++Q1ngUuSNryolUruOgimDcPbropX+8qIlK4zAt4akx5eblnc1N4d9hxR6iqioXi1loray8tIlIwzGyGu5fXd1yLvRK4Lmbwhz9EAFx3XdrViIikq6QCAGD33WGPPSIIPv447WpERNJTcgEAMRZQXQ1XXZV2JSIi6SnJABgwACoq4NJLYXHehqFFRApLSQYAwO9/H11Al16adiUiIuko2QD4/vfh//2/6AZ65520qxERyb+SDQCA886Li8IuuijtSkRE8q+kA2DTTeHww+H662P/YBGRUlLSAQCxZWSrVrFWkIhIKSn5AOjRA44+GsaNgzlz0q5GRCR/Sj4AAM44A9ZZJxaKExEpFQoAYt/gE0+EO++E559PuxoRkfxQACROPhk6d9YG8iJSOhQAifXXj66gSZPg3/9OuxoRkdxTANRyzDGw8cZw1lnaQF5EWj4FQC3rrBNdQP/6Fzz6aNrViIjklgJgJb/8JfTuDWefrbMAEWnZFAAradculoiYMQPuvTftakREckcBUIdDDoF+/aI7aNmytKsREckNBUAdWreOpSFmz4Zbbkm7GhGR3FAArMYPfwjbbhvdQV99lXY1IiLZpwBYjVatYpno+fPhhhvSrkZEJPsUAGuw996wyy6xe9hnn6VdjYhIdtUbAGa2lpk9Y2YvmtksMzs/ae9jZtPNbK6Z3W5m7ZL29snX85Lv9671Wmcm7a+a2d65+qGyxSzOAhYuhGuvTbsaEZHsasgZwJfAYHffGugPDDWzgcAo4Ap37wssAY5Ijj8CWOLumwJXJMdhZlsABwHfA4YCfzaz1tn8YXJhl11g6FC4+GL46KO0qxERyZ56A8DDJ8mXbZObA4OBu5L2ccD+yeOK5GuS7w8xM0vax7v7l+4+H5gHDMjKT5FjF14IixfD5ZenXYmISPY0aAzAzFqb2QvAImAy8F/gA3dfmhxSBXRPHncH3gJIvv8hsEHt9jqeU9C22w5+/GMYPRreey/takREsqNBAeDuy9y9P9CD+Ku9X12HJfe2mu+trv0bzGykmVWaWWV1dXVDysuLCy6IgeCLL067EhGR7GjULCB3/wB4HBgIdDKzNsm3egALksdVQE+A5PvrA4trt9fxnNrvMdbdy929vKysrDHl5VS/fnDooTEY/PbbaVcjItJ8DZkFVGZmnZLHawN7ALOBx4AfJ4eNAO5PHk9Ivib5/lR396T9oGSWUB+gL/BMtn6QfDj33Fga4sIL065E0vbb38IDD6RdhUjzNOQMoBvwmJnNBJ4FJrv7g8DpwElmNo/o46+5XOoGYIOk/STgDAB3nwXcAbwCPAwc7e5FtdJOnz4wciT89a/w3/+mXY2k5ZlnYnrwRRelXYlI85gX8JrH5eXlXllZmXYZ37BwIWyyCfzoR3DzzWlXI2n40Y/gnnvi8YIF0K1buvWIrMzMZrh7eX3H6UrgRurWDY49NhaJmzUr7Wok3157LZYJz2Ti6wcfTLcekeZQADTBaafBuuvCOeekXYnk22WXxZ4R118fXYITJqRdkUjTKQCaYIMN4JRT4i/BZ59NuxrJl4ULYdw4+PnPY+/oTAb++U/49NO0KxNpGgVAE51wAmy4YWwdKaXhqqtg6dIIf4gA+OILmDw53bpEmkoB0ETrrgtnnhn/+B9/PO1qJNc+/BDGjIkB4E03jbZddoFOndQNJMVLAdAMRx0F3btrA/lScP31sRjg6aevaGvbFvbdNwaCtXWoFCMFQDOsvTb87nfw1FMwaVLa1UiufPklXHklDBkS60LVlslAdTVMm5ZObSLNoQBopl/8Iq4LOPtsWL487WokF26+OQaAa//1X2Po0DgTUDeQFCMFQDO1bQvnnw8vvgh33pl2NZJty5bBpZfCNtvAHnus+v3114fddlMASHFSAGTBQQfBlltGd9DSpfUfL8Xj/vvj4q/TT48d4uqSycCcOXGcSDFRAGRB69axQNxrr8Hf/552NZIt7jBqFHznOzH7Z3VqrgrWWYAUGwVAlmQyMGAAnHdeDBpK8XviiVj47ZRToE2b1R/Xqxf0768AkOKjAMiSmg3k33orpgxK8bvkEigriyt/65PJwH/+ox3jpLgoALJoyBDYffcIAi0PUNxmzoSHHoLjjovpvvWpqIhZYBMn5r42kWxRAGRRzVnAokVw9dVpVyPNcckl0KED/OY3DTt+m23iokB1A0kxUQBk2aBBsN9+8QGyZEna1UhTvPEGjB8fm/906dKw55hFN9Ajj8T6QCLFQAGQA7//PXzwQSwdLMXn8svjA/3EExv3vEwmuv6mTs1NXSLZpgDIga23jmsDrroK3n037WqkMd57L7b8/NnPoGfPxj13992hY0d1A0nxUADkyPnnR1fAH/+YdiXSGNdcA59/Hpv+NFb79rE0xIQJWhZEioMCIEc22yzWCRozBt58M+1qpCE+/RT+9KcYw9lii6a9RiYT6wbNmJHd2kRyQQGQQzVbRl5wQbp1SMPccAMsXlz3om8NNWxYXBmubiApBgqAHOrVK/YMuOkmrRNT6L7+GkaPhp12iltTdekCO++sAJDioADIsTPPhLXWgnPPTbsSWZPbb4+uuub89V8jk4kLyebPb/5rieRSvQFgZj3N7DEzm21ms8zs+KT9PDN728xeSG771nrOmWY2z8xeNbO9a7UPTdrmmdkZufmRCkvXrrF/8PjxsWS0FB73uG5jiy2iC6e5Kiri/oEHmv9aIrnUkDOApcDJ7t4PGAgcbWY1Q2RXuHv/5DYJIPneQcD3gKHAn82stZm1Bq4F9gG2AA6u9Tot2imnxN6xv/1t2pVIXR5+GF56KWb+tMrCOfEmm0SYqBtICl29v+7uvtDdn0sefwzMBrqv4SkVwHh3/9Ld5wPzgAHJbZ67v+7uXwHjk2NbvE6d4sPlwQfh6afTrkZWNmoU9OgBBx+cvdfMZGI10Q8+yN5rimRbo/7eMbPewDbA9KTpGDObaWY3mlnnpK078Fatp1UlbatrX/k9RppZpZlVVldXN6a8gnbccbDRRnDWWdpAvpBMnx4f1CedBO3aZe91M5nYHOihh7L3miLZ1uAAMLOOwN3ACe7+ETAG2AToDywERtccWsfTfQ3t32xwH+vu5e5eXlZW1tDyCl6HDtEF9PjjcOONCoFCMWpUnKEdeWR2X3eHHSLw1Q0khaxBAWBmbYkP/1vc/R4Ad3/X3Ze5+3LgL0QXD8Rf9rUvou8BLFhDe8kYOTI+GI48EoYPj0XHJD2vvgr33QdHHw3rrpvd127VKi4omzQJvvoqu68tki0NmQVkwA3AbHe/vFZ7t1qHHQC8nDyeABxkZu3NrA/QF3gGeBboa2Z9zKwdMVBcUn8ftW8P//53LDb2xBPwve/FgnHaRzgdl14a/0+OOy43r5/JwEcfwZNP5ub1RZqrIWcAOwGHAoNXmvJ5iZm9ZGYzgd2BEwHcfRZwB/AK8DBwdHKmsBQ4BniEGEi+Izm2pLRpE6tMvvJKbCBz6qlQXh5bD0r+LFgAN98cy3VstFFu3mOPPWIzGXUDSaEyL+DO6PLycq+srEy7jJxxh3vvhWOPjfVjjj46NpRZb720K2v5TjstrvydOzc2fc+Vigp44YXo7rO6RsFEcsDMZrh7eX3H6UrgFJnBD38Is2fDMcfAtddCv35w990aJM6lDz6A666Dn/wktx/+EN1Ab74ZVwaLFBoFQAFYb73YQnLatNiE/Mc/jr8ctYpoblx3HXz8cXaWfajP8OER9OoGkkKkACggAwZAZWUMDE+ZEleTXnGFBomz6Ysv4MorYc89Yx/fXOvaFQYOVABIYVIAFJg2beDkk2HWLNh117hAqSYYpPluvjl2acvHX/81Mpn4//f22/l7T5GGUAAUqN69Y+mIO++Ed96J6weOPz66LqRpli2LqZ/bbQeDB+fvfbU4nBQqBUABM4vxgNmz4de/jt2q+vWLi5ek8e67L2b9nH56fmfkfPe7sOmm6gaSwqMAKALrrx8zhJ56KjYcOeAA2H9/eOut+p8rwT2Wfdh005h5lU9m0Q00ZYrO4KSwKACKyMCBsdfsqFHw6KMxSHzVVdG1IWv2+OPw7LOxNHfr1vl//0wmloR49NH8v7fI6igAikzbtnER06xZsfXgCSfE+MBzz6VdWWEbNSqu+D3ssHTef6ed4uxN3UBSSBQARapPn1hobPx4qKqC7bePGUOffJJ2ZYXnhRfgkUdiEH3ttdOpoU2b2G3swQc1rVcKhwKgiJnBT38Kc+bESqNXXKGdqOpyySXQsSMcdVS6dWQysHhxjOWIFAIFQAvQqROMGQP/+U9cVVxRAT/6keadQ2zMfvvt8KtfQefO9R+fS3vvHZvOKKClUCgAWpAdd4yxgD/+MbqH+vWDa64p7UHi0aNj0PfEE9OuJPYcGDwY7r9faz1JYVAAtDDt2sEZZ8DLL8OgQbHS6KBB0Q9eaqqrY/e1Qw6B7mvaxTqPMhmYNy+67UTSpgBooTbZBB5+GG69Ff73v9hz4JRT4NNP064sf/70J/j889hzoVDst1/cqxtICoECoAUzg4MPjr82Dz88ukO22AImTky7stz75JPo/qqoiK6wQtGjRyxFoQCQQqAAKAGdO8PYsfCvf8VsmOHDYy38BS14R+YbboAlS/K76FtDZTLw9NOxKJ1ImhQAJWTnneH55+HCC2Nhsn794M9/bnmDxF9/Hfsu77JLjH8UmoqKGAQuhTMxKWwKgBLTrh2cfXYMEm+/fWxDudNOLWvHqvHjYzOdQvzrH2CrraBXL3UDSfoUACVq001h8uRYH/+//4Vtt4Xzziv+swH3uPBryy1h333TrqZuNYvDPfoofPZZ2tVIKVMAlDCzmCI5Z04MFp9/fsxTL+YLyCZNirOb004r7E3YM5mYoTRlStqVSClTAAgbbBBnAuPGxWqjW29dvP3To0ZBz55w0EFpV7Jmu+4aV22rG0jSVG8AmFlPM3vMzGab2SwzOz5p72Jmk81sbnLfOWk3M7vazOaZ2Uwz27bWa41Ijp9rZiNy92NJUxx2WARAjx4xU+ikk2IJ42Lx9NMx0+mkk2LV1ELWrh3ss08Mxi9fnnY1UqoacgawFDjZ3fsBA4GjzWwL4Axgirv3BaYkXwPsA/RNbiOBMRCBAZwL7AAMAM6tCQ0pHJtvDtOmwTHHxOJyO+0UYwTFYNSomPJ65JFpV9IwmUxMBX3mmbQrkVJVbwC4+0J3fy55/DEwG+gOVADjksPGAfsnjyuAv3uYBnQys27A3sBkd1/s7kuAycDQrP40khVrrRVX0d5zTyxbsM02MbOmkM2eHWvsHHNMXOtQDPbZJ5aJVjeQpKVRYwBm1hvYBpgOdHX3hRAhAWyUHNYdqL1ZYVXStrp2KVAHHBBrCH3/+zFIfOSRhTtr5dJLY63/Y49Nu5KG69wZfvCDCC6RNDQ4AMysI3A3cIK7f7SmQ+to8zW0r/w+I82s0swqq6urG1qe5Mi3vx3bKZ51Viystv32McumkFRVwT/+EctdlJWlXU3jZDLwyitxpiWSbw0KADNrS3z43+Lu9yTN7yZdOyT3i5L2KqBnraf3ABasof0b3H2su5e7e3lZsf1rbqHatoWLLopdtd5/P0Jg7NjCWdL4qqtiIPXkk9OupPEymbh/4IF065DS1JBZQAbcAMx298trfWsCUDOTZwRwf632w5LZQAOBD5MuokeAvcysczL4u1fSJkVizz3hxRdjiYVf/SqmWn74Ybo1ffABXH89HHhgbJNZbPr0iS42jQNIGhpyBrATcCgw2MxeSG77AhcDe5rZXGDP5GuAScDrwDzgL8BvANx9MfB74NnkdkHSJkWka9dYZvqPf4S7744B4jRnsYwZAx9/HBd+FatMJqavLta/Bskz80I5j69DeXm5V1ZWpl2GrMbTT8fg8NtvRyCcdBK0yuOlhV98Ab17Q//+EUrF6plnYIcd4mK8Qw5JuxppCcxshruX13ecrgSWJhs0KFYXzWRi05Xhw2MXrnwZNy7m0Rfqom8NVV4O3bqpG0jyTwEgzdK5M9x1F1x7LUydGstIPPZY7t932TK47LIYkN5tt9y/Xy61ahU7hT30EHz5ZdrVSClRAEizmcFvfgPTp8f6NkOGwLnnwtKluXvPmovUTj+9sBd9a6hMJnYxe/zxtCuRUqIAkKzZemuorIQRI+CCC2Jl0aqq7L+Peyz70Lcv7L9//ccXgyFDYJ111A0k+aUAkKzq2BH+9rcY0HzuuQiFbM9xnzo1Fq079VRo3Tq7r52WtdaCvfeOACjgeRnSwigAJCcOOSQCoFev6N448cTs9W+PGhXTUQ89NDuvVygymThjev75tCuRUqEAkJzZbLNYWfS44+DKK2HHHZu/5MFzz8VOZiecEH81tyTDhsWAsLqBJF8UAJJT7dvHUg333Qfz58eFY7fe2vTXu+QSWHdd+PWvs1djoSgri5DU4nCSLwoAyYuKilhZdOut4Wc/gyOOgE8/bdxrvP463HlnfPh36pSbOtOWycR/pzffTLsSKQUKAMmbXr1imuPZZ8dA8fbbw0svNfz5o0fH+vknnJCzElNXURH3WhxO8kEBIHnVpg1ceGH04y9ZAgMGxGJu9c18WbQolqM+9FD41rfyU2saNtssdmXTOIDkgwJAUjFkSHR17LprdOn89Kexsufq/OlPMYvo1FPzV2NaMpm4mjrtlVal5VMASGq6doVJk2Ja5733xgDx9OmrHvfJJ7HUxP77x1/HLV0mA19/HfsviOSSAkBS1apVLOX8r3/F1zvvHNs7Ll++4pi//CW6i4p90beGGjQINtxQ3UCSewoAKQgDB8YFUBUVEQjDhkW//1dfweWXR1fRDjukXWV+tG4dK6tOnBhnAiK5ogCQgtGpU0zzHDMm+sC33jq2eayqKp2//mtkMjEm8u9/p12JtGQKACkoZjEo/MwzEQjXXANbbQVDh6ZdWX7ttVdcRKduIMklBYAUpK22ipVFzzknpom2hCWfG6NDB9hjj7gqWIvDSa4oAKRgdegQy0oPHJh2JenIZGL5jFmz0q5EWioFgEiB2m+/uFc3kOSKAkCkQHXrFldKa3E4yRUFgEgBy2RiQHzhwrQrkZZIASBSwDKZuH/wwXTrkJap3gAwsxvNbJGZvVyr7Twze9vMXkhu+9b63plmNs/MXjWzvWu1D03a5pnZGdn/UURani23hD591A0kudGQM4CbgLpmYV/h7v2T2yQAM9sCOAj4XvKcP5tZazNrDVwL7ANsARycHCsia2AWZwH//Gfj908QqU+9AeDuTwKLG/h6FcB4d//S3ecD84AByW2eu7/u7l8B45NjRaQemUyshDp5ctqVSEvTnDGAY8xsZtJF1Dlp6w68VeuYqqRtde2rMLORZlZpZpXV1dXNKE+kZdhll7gqWtNBJduaGgBjgE2A/sBCYHTSXtf1mr6G9lUb3ce6e7m7l5eVlTWxPJGWo21b2HffGAhetiztaqQlaVIAuPu77r7M3ZcDfyG6eCD+su9Z69AewII1tItIA2QyUF0N06alXYm0JE0KADPrVuvLA4CaGUITgIPMrL2Z9QH6As8AzwJ9zayPmbUjBop1QivSQEOHxpmAuoEkmxoyDfQ24GlgczOrMrMjgEvM7CUzmwnsDpwI4O6zgDuAV4CHgaOTM4WlwDHAI8Bs4I7kWBFpgPXXh91203RQyS7zAl5qsLy83CsrK9MuQ6QgXHMNHHsszJlTGltjStOZ2Qx3L6/vOF0JLFIkaq4KfuCBdOuQlkMBIFIkevWC/v3VDSTZowAQKSKZDDz1VMwIEmkuBYBIEclkYPlymDQp7UqkJVAAiBSRbbeF7t3VDSTZoQAQKSI1i8M98gh88UXa1UixUwCIFJlMBj77DKZOTbsSKXYKAJEis/vu0LGjrgqW5lMAiBSZ9u1jaYgJE2JAWKSpFAAiRSiTiX2CZ8xIuxIpZgoAkSI0bBi0bq1uIGkeBYBIEerSBXbeWdNBpXkUACJFKpOBl16C+fPTrkSKlQJApEhVJLtqa3E4aSoFgEiR2mQT2GILdQNJ0ykARIpYJgNPPAFLlqRdiRQjBYBIEctkYqP4hx9OuxIpRgoAkSK2ww6w0UbqBpKmUQCIFLFWrWC//eChh+Crr9KuRoqNAkCkyGUy8NFH8OSTaVcixUYBIFLk9tgD1l5bVwVL4ykARIrcOuvAnnvGOIB72tVIMak3AMzsRjNbZGYv12rrYmaTzWxuct85aTczu9rM5pnZTDPbttZzRiTHzzWzEbn5cURKUyYDb74JM2emXYkUk4acAdwEDF2p7Qxgirv3BaYkXwPsA/RNbiOBMRCBAZwL7AAMAM6tCQ0Rab7hw2O3MHUDSWPUGwDu/iSweKXmCmBc8ngcsH+t9r97mAZ0MrNuwN7AZHdf7O5LgMmsGioi0kRdu8LAgZoOKo3T1DGAru6+ECC53yhp7w68Veu4qqRtde0ikiWZTOwPUFWVdiVSLLI9CGx1tPka2ld9AbORZlZpZpXV1dVZLU6kJctk4v7mm7VTmDRMUwPg3aRrh+R+UdJeBfSsdVwPYMEa2lfh7mPdvdzdy8vKyppYnkjp6dcPttkGzjoLNt4Yfv5zuPNO+PDDtCuTQtXUAJgA1MzkGQHcX6v9sGQ20EDgw6SL6BFgLzPrnAz+7pW0iUiWmMFjj8Gtt8a00AkT4MADYcMNYfBgGD0a5szRVFFZwbye3wYzuw3YDdgQeJeYzXMfcAfQC3gT+Im7LzYzA64hBng/A37h7pXJ6xwOnJW87EXu/rf6iisvL/fKysom/FgisnQpTJsGEyfCgw/Cy8lE7k02iS0lhw+HH/wgNpmXlsXMZrh7eb3H1RcAaVIAiGTP//4XYTBxIkydCl98AR06xNnC8OGw777QrVvaVUo2KABEZLU++yxCoObsoGbm0HbbxdnBsGFQXh6LzUnxUQCISIO4x97CNWEwbVrMIuraFfbZJ84O9twT1lsv7UqloRQAItIk770XG8xMnBj3H3wAbdvCLrtEGAwbBpttlnaVsiYKABFptqVL4amnVpwdvPJKtPftu6Kr6Ac/gHbt0q1TvkkBICJZN3/+ioHkxx6DL7+EddddMZC8zz5xDYKkSwEgIjn16acwZUqcGUycCAuSSzvLy1d0FW27rQaS06AAEJG8cYcXX1wRBtOnR9vGG8f00mHDYK+9oGPHtCstDQoAEUlNdXXsU1wzkPzRR7FxzQEHwCGHxC5mbdqkXWXL1dAA0MmZiGRdWRkcdhjcfnvMKpo6NT74J06McYLu3eGEE6CyUktTpEkBICI51bYt7L47XH89vPMO3HMP7LwzjBkD228fi9hdeGEMMEt+KQBEJG/at49uoLvvjjAYOzYuODvnHPjOdyIYrrsO3n8/7UpLgwJARFLRuTP88pfwxBPwxhvwhz/A4sVw1FGxJtH++8Ndd8WaRZIbCgARSd23vw1nngmzZsFzz8Gxx8ZMop/8JGYSHXkkPP64NrrJNgWAiBQMs9jUZvToWKDu0UehogLGj49xhN694YwzVixtLc2jABCRgtS6dVxhPG4cvPtubHTz/e/DZZfFff/+8fjtt9OutHgpAESk4HXoAAcfvOKK46uvjgHlU0+Fnj3juoKbborrDaThFAAiUlQ22mjFGMGrr8YMovnz4Re/iBlFBx0UVyR//XXalRY+BYCIFK3NNoPzz4d582LV0sMPh3/+E/bbD771LTjmmNjfQBeb1U0BICJFzwwGDYJrr40uogkTYPBguOGGaK8dFLKCAkBEWpR27eIM4Pbb42KzG2+EXr0iAPr2XREU1dVpV5o+LQaVqkxcAAAG/klEQVQnIiWhqgpuuw3+8Q+YOTMWoxs6NNYm2nzzOEvo3r1lLF+t1UBFRFZj5ky45Za41Z5GuvbacZaw2War3jbYIL16G0sBICJSD/cYM3jttVVvr78eW2LW6NKl7mDYdNOYplpI8hIAZvYG8DGwDFjq7uVm1gW4HegNvAEc6O5LzMyAq4B9gc+An7v7c2t6fQWAiKTl669jjaK6wqGq6pvH9ugRYbDy2UOfPrEaar41NACysSXD7u7+Xq2vzwCmuPvFZnZG8vXpwD5A3+S2AzAmuRcRKTht28YHet++saNZbZ9+GjOKVg6GO+6AJUtWHNe6daxyWteZQ/fuMXspTbnYk6cC2C15PA54nAiACuDvHqcc08ysk5l1c/eFOahBRCRnOnSArbeO28ref7/us4apU+Hzz1cct846qx9v6NIlPz9HcwPAgUfNzIHr3X0s0LXmQ93dF5rZRsmx3YG3aj23KmlTAIhIi7HBBjHVdNCgb7YvXx4DzisHw/PPxyY5y5Z98zX23DNmLeVScwNgJ3dfkHzITzazOWs4tq6TnVUGIMxsJDASoFevXs0sT0SkMLRqFesW9ewJQ4Z883tffRXLWcyduyIY8nEW0KwAcPcFyf0iM7sXGAC8W9O1Y2bdgEXJ4VVAz1pP7wEsqOM1xwJjIQaBm1OfiEgxaNcurkXYfPP8vm+TL3kwsw5mtm7NY2Av4GVgAjAiOWwEcH/yeAJwmIWBwIfq/xcRSU9zzgC6AvfG7E7aALe6+8Nm9ixwh5kdAbwJ/CQ5fhIxBXQeMQ30F814bxERaaYmB4C7vw6sMgbu7u8DQ+pod+Dopr6fiIhkVwtY9UJERJpCASAiUqIUACIiJUoBICJSohQAIiIlqqCXgzazauB/zXiJDYH36j0q/1RX46iuxlFdjdMS6/q2u5fVd1BBB0BzmVllQ5ZEzTfV1Tiqq3FUV+OUcl3qAhIRKVEKABGREtXSA2Bs2gWshupqHNXVOKqrcUq2rhY9BiAiIqvX0s8ARERkNVpcAJhZTzN7zMxmm9ksMzs+7ZpqM7PWZva8mT2Ydi01ku057zKzOcl/t0H1Pyv3zOzE5P/hy2Z2m5mtlWItN5rZIjN7uVZbFzObbGZzk/vOBVLXpcn/y5lmdq+ZdSqEump97xQzczPbsFDqMrNjzezV5PftkkKoy8z6m9k0M3vBzCrNbEC237fFBQCwFDjZ3fsBA4GjzWyLlGuq7XhgdtpFrOQq4GF3/y6xwmvq9ZlZd+A4oNzdtwRaAwelWNJNwNCV2s4Aprh7X2BK8nW+3cSqdU0GtnT3rYDXgDPzXRR114WZ9QT2JJaKT8NNrFSXme1O7Fm+lbt/D7isEOoCLgHOd/f+wO+Sr7OqxQWAuy909+eSxx8TH2bd060qmFkPYBjw17RrqWFm6wE/AG4AcPev3P2DdKv6P22Atc2sDbAOdewgly/u/iSweKXmCmBc8ngcsH9ei6Luutz9UXdfmnw5jdh9L/W6ElcAp1HHdrD5sJq6jgIudvcvk2MWrfLEdOpyYL3k8frk4Pe/xQVAbWbWG9gGmJ5uJf/nSuKXf3nahdTyHaAa+FvSNfXXZIe3VLn728RfYm8CC4kd5B5Nt6pVdK3Z1S653yjleupyOPBQ2kUAmFkGeNvdX0y7lpVsBuxiZtPN7Akz2z7tghInAJea2VvEv4Wsn8m12AAws47A3cAJ7v5RAdQzHFjk7jPSrmUlbYBtgTHuvg3wKel0ZXxD0p9eAfQBvgV0MLND0q2quJjZ2USX6C0FUMs6wNlEV0ahaQN0JrqMTyV2NLR0SwLizOREd+8JnEhylp5NLTIAzKwt8eF/i7vfk3Y9iZ2AjJm9AYwHBpvZP9ItCYAqoMrda86S7iICIW17APPdvdrdvwbuAXZMuaaVvWtm3QCS+7x3HayOmY0AhgM/88KY670JEeYvJv8GegDPmdnGqVYVqoB7PDxDnKHnfYC6DiOI33uAOwENAtcnSe4bgNnufnna9dRw9zPdvYe79yYGM6e6e+p/0br7O8BbZrZ50jQEeCXFkmq8CQw0s3WS/6dDKIDB6ZVMIP6Rktzfn2It/8fMhgKnAxl3/yztegDc/SV338jdeyf/BqqAbZPfv7TdBwwGMLPNgHYUxuJwC4Bdk8eDgblZfwd3b1E3YGdi8GQm8EJy2zftulaqcTfgwbTrqFVPf6Ay+W92H9A57ZqSus4H5gAvAzcD7VOs5TZiLOJr4sPrCGADYvbP3OS+S4HUNQ94q9bv/3WFUNdK338D2LAQ6iI+8P+R/J49BwwukLp2BmYALxLjmNtl+311JbCISIlqcV1AIiLSMAoAEZESpQAQESlRCgARkRKlABARKVEKABGREqUAEBEpUQoAEZES9f8BasBmrOL45YAAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "#再用图看一下\n",
    "plt.plot(Ks2, np.array(chis), 'b-')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x12555710>]"
      ]
     },
     "execution_count": 32,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD8CAYAAACMwORRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAHaRJREFUeJzt3XuY1VXd/vH3h2OaKRhYCijnk8TJkTTTBCUBc0YfMSEzTfBAgolHeEJE8DEPyMGkzMOvrFRSMB3AQstjqMSgqBwEERFQDNQ8pCkin98fa9MM48DsYfaetfd336/rmgv2nq+z70uGmzXru/Za5u6IiEiy1IsdQEREMk/lLiKSQCp3EZEEUrmLiCSQyl1EJIFU7iIiCaRyFxFJIJW7iEgCqdxFRBKoQawXbtasmbdu3TrWy4uI5KXFixe/7e7Nq7suWrm3bt2asrKyWC8vIpKXzOz1dK7TtIyISAKp3EVEEkjlLiKSQCp3EZEEUrmLiCSQyl1EJIFU7iIiCZR35f788zB2LOh0QBGRncu7cl+wAK69Fv7619hJRERyV96V+9lnQ6tWcMUVGr2LiOxM3pV748YwfjwsXAjz5sVOIyKSm/Ku3AHOOAPatQuj923bYqcREck9eVnuDRvClVfCkiVw//2x04iI5J68LHeAH/wAOncOUzSffx47jYhIbsnbcq9fHyZOhBUr4J57YqcREckteVvuACefDD16wIQJ8NlnsdOIiOSOvC73evVg0iR49VW4887YaUREckdelzvA974HffqEKZpPP42dRkQkN+R9uZvB1VfD+vVw222x04iI5Ia8L3eAY4+Fo46C//s/+Pjj2GlEROJLRLmbhbn3t96CX/4ydhoRkfgSUe4QRu7f/W7YVOzDD2OnERGJKzHlDmH0/s47MH167CQiInElqtz79IHiYpg8Gf71r9hpRETiSVS5Q1gS+f77MGVK7CQiIvEkrtx79IDvfx+mTYPNm2OnERGJI61yN7MBZrbSzFab2ZgqPj/VzJakPlaZ2XuZj5q+CRPCksjrr4+ZQkQknmrL3czqAzOAgUBXYKiZda14jbuPdvee7t4T+AUQdSPeLl3gtNPg5pth48aYSURE4khn5N4HWO3ua9x9CzATKNnF9UOB6Ps0Xnll2EzsmmtiJxERqXvplHsLYH2FxxtSz32BmR0EtAEerX202mnXDs46C379a3j99dhpRETqVjrlblU8t7OjqYcAs9y9yuMzzOwcMyszs7LNdXC3c9y48r1nREQKSTrlvgFoVeFxS+DNnVw7hF1Mybj7re5e5O5FzZs3Tz/lbjrwQDj3XPjNb2D16qy/nIhIzkin3BcBHcysjZk1IhR4aeWLzKwT0BR4JrMRa2fsWGjUCK66KnYSEZG6U225u/tWYCQwH1gB3Ovuy8xsopkVV7h0KDDT3Xc2ZRPF/vvDyJFw112wfHnsNCIidcNidXFRUZGXlZXVyWu9/Ta0aQMDBsB999XJS4qIZIWZLXb3ouquS9w7VKvSrBmMHg2zZsHzz8dOIyKSfQVR7gAXXQRNmsD48bGTiIhkX8GUe5MmcOmlMHcuLFwYO42ISHYVTLkDXHBBmKK54orYSUREsqugyn2vvcLSyEcegSeeiJ1GRCR7CqrcAUaMCMsjr7gCcmvRpohI5hRcue+xR9iW4KmnwgheRCSJCq7cAYYNC1sTjBun0buIJFNBlnvjxmFJ5KJFMGdO7DQiIplXkOUO8KMfQfv2Ye5927bYaUREMqtgy71hw3Ac34svhneuiogkScGWO8CQIdC1a5ii2bo1dhoRkcwp6HKvXx8mToSVK+Huu2OnERHJnIIud4CTToJevcIUzWefxU4jIpIZBV/u9erBpEnw2mvhxCYRkSQo+HIHGDQIDjsslPwnn8ROIyJSeyp3yg/R3rABbr01dhoRkdpTuaf06wdHHw3XXAMffxw7jYhI7ajcU8zCtMw//wkzZsROIyJSOyr3Cr797XDO6nXXwQcfxE4jIrL7VO6VTJoE77wD06fHTiIisvtU7pUUFcGJJ8LkyfDuu7HTiIjsHpV7Fa66KkzL3Hhj7CQiIrtH5V6F7t3h1FPD1MymTbHTiIjUnMp9JyZMgP/8J9xcFRHJNyr3nejcGU4/HX75S3jjjdhpRERqRuW+C9u3Ar7mmthJRERqRuW+C23bhvNWb7sN1q6NnUZEJH0q92qMG1e+c6SISL5Iq9zNbICZrTSz1WY2ZifXfN/MlpvZMjNLzNEXLVvCeefBnXfCqlWx04iIpKfacjez+sAMYCDQFRhqZl0rXdMBGAsc4e4HAxdmIWs0Y8ZA48Zh/buISD5IZ+TeB1jt7mvcfQswEyipdM3ZwAx3/xeAuydqdfjXvw6jRsE998CyZbHTiIhUL51ybwGsr/B4Q+q5ijoCHc1sgZk9a2YDMhUwV1x6Key1F1x5ZewkIiLVS6fcrYrnvNLjBkAH4GhgKHC7mTX5whcyO8fMysysbPPmzTXNGtVXvwoXXQSzZ8Nzz8VOIyKya+mU+wagVYXHLYE3q7jmQXf/zN1fA1YSyn4H7n6ruxe5e1Hz5s13N3M0o0dD06Zh/buISC5Lp9wXAR3MrI2ZNQKGAKWVrnkA6AtgZs0I0zRrMhk0F+yzD1x2GcybB888EzuNiMjOVVvu7r4VGAnMB1YA97r7MjObaGbFqcvmA++Y2XLgMeBSd38nW6FjGjUK9tsPrrgidhIRkZ0z98rT53WjqKjIy8rKorx2bU2dGubfH30U+vaNnUZEComZLXb3ouqu0ztUd8N558EBB4TRe6R/G0VEdknlvhv22CNsS7BgAcyfHzuNiMgXqdx307BhcNBBoeQ1eheRXKNy302NGoU3NC1eDA8+GDuNiMiOVO61cPrp0KFDmHvfti12GhGRcir3WmjQIGwmtnQp3Htv7DQiIuVU7rV06qnQrVuYotm6NXYaEZFA5V5L9erBxIlhr/e77oqdRkQkULlnwIknQu/eYfS+cmXsNCIiKveMMINp0+Bf/wpTNJdcAu+/HzuViBQylXuGHHlkmJo580yYMiWsorn9dvj889jJRKQQqdwz6Gtfg9tug7Iy6NQJzj4bDj0UnnoqdjIRKTQq9yzo3RuefBJmzoS334ajjoIhQ2DdutjJRKRQqNyzxCwsk3z5ZZgwAUpLw2h+wgT4+OPY6UQk6VTuWbbnnmEVzcsvQ0lJeNNT585hVK89aUQkW1TudeTAA0OhP/kkNGsGQ4eGm7CLF8dOJiJJpHKvY0ceCYsWhRuvr7wSbrgOHw7//GfsZCKSJCr3COrXD4W+alU40enOO6FjR5g8GbZsiZ1ORJJA5R7RPvuEQl+6NIzoL700vAlq7lzNx4tI7ajcc0CnTqHQH3oo7FVzwgkwcCCsWBE7mYjkK5V7Dhk4EF56KRzA/eyz0L07jB4dtjUQEakJlXuOadgQLrww3GwdNgymTw/z8b/+tbYyEJH0qdxzVPPmcMst8NxzcPDBcN554Z2vjz8eO5mI5AOVe47r2RMeewzuuy/sNNm3L5xyCqxdGzuZiOQylXseMIPBg8MN1kmTwo3Xzp3D2a0ffRQ7nYjkIpV7HtljDxg3LhwIMngwXH11WGlz111aOikiO1K556GWLeEPf4AFC2D//eGHP4QjjgjvfBURAZV7XvvWt2DhQvjNb2DNGujTB378Y9i4MXYyEYlN5Z7n6tULpz+tWgWXXw533x2WTl53HXz6aex0IhJLWuVuZgPMbKWZrTazMVV8/kwz22xmS1IfwzMfVXZl773h2mth2TI45hgYMyYsoXzwQc3HixSiasvdzOoDM4CBQFdgqJl1reLSP7p7z9TH7RnOKWlq3x4eeAAefhgaN4YTTwxlf/fd8MEHsdOJSF1JZ+TeB1jt7mvcfQswEyjJbiyprf79YckSuOmmMGVz2mnhjVHFxfC738F778VOKCLZlE65twDWV3i8IfVcZSeb2YtmNsvMWmUkndRKw4YwalQ4u3XBAjj//FD4Z5wB++0Hxx8fbsa++27spCKSaemUu1XxXOVZ3DlAa3fvDvwVuLPKL2R2jpmVmVnZ5s2ba5ZUdlu9emFlzZQp8PrrYYXNhRfC8uVw1lnwta/BccfB7beHA71FJP+ZV3O3zcwOBya4+3Gpx2MB3P3nO7m+PvCuu++zq69bVFTkZWVluxVaMsM97F0za1bY3uDVV8NBIkcfHd4kddJJofhFJHeY2WJ3L6ruunRG7ouADmbWxswaAUOA0kovtn+Fh8WAdiLPA2ZwyCHw85+HXSiffz6sslm/HkaMgAMOCHvZzJihtfMi+abacnf3rcBIYD6htO9192VmNtHMilOXXWBmy8zsBeAC4MxsBZbsMAublF19Nbz8cthXftw42LQJRo6EFi3CaVHTp8OGDbHTikh1qp2WyRZNy+SP5cth9uwwffPii+G5ww8PUzcnnwwHHRQ3n0ghSXdaRuUuNbJqVSj5WbPCNA7AoYeGbYhPPhnato2bTyTpVO6Sda++Gkb0990H2/8oe/cOI/rBg6FDh7j5RJJI5S51au3a8qmbZ58Nz/XoUV70nTtHjSeSGCp3iWbdOrj//lD0CxaE5w4+OEzdDB4MXbuGG7giUnMqd8kJb7xRXvRPPRXW1nfuXD6i795dRS9SEyp3yTlvvQV/+lMo+scfh23boHXrsN9NSUlYatmwYeyUIrlN5S45bdMmKC0NWxL/9a/wySfQpAkMGhTKfuDAsI2xiOxI5S5546OP4JFHQtnPmRP2t2nYMGyDUFISyr6VtqITAVTukqc+/xyeeaZ8VL9qVXi+V6/you/ZU/P0UrhU7pIIK1eGki8thaefDjdkW7Uqn6f/znegUaPYKUXqjspdEmfTJpg3L5T9ww/Df/4T5uUHDgxFP3BgmLcXSTKVuyTaxx/D3/4Win7OnFD8DRqEkXxxcfho3Tp2SpHMU7lLwdi2LRxAsn2efkVqw+nu3cvn6Q85RPP0kgwqdylYr7wSir60FP7+91D+LVrACSeEsu/bNxweLpKPVO4ihGWV8+aFop8/Pyy73GsvGDAgFP2gQbDvvrFTiqRP5S5SySefwKOPls/Tb9wYjhU88sjy1TfaslhyncpdZBe2bQvbFG+fp1+6NDzfrRuceCL85Cew//67/hoiMWTyDFWRxKlXD/r0CccKvvRS2Jt+6lRo3jycKdu2LVx8cViFI5KPVO4ihDK/8MIwbbNyJZx6KkybFp4fOxbeeSd2QpGaUbmLVNKuHfz2t2FJZUkJXHcdtGkD48fDe+/FTieSHpW7yE507Ah33RWmbY47DiZNCm+MmjQJPvggdjqRXVO5i1Tj4IPDObFLloSdKsePDyP5a6+Ff/87djqRqqncRdLUowc88EBYZXPYYWEuvm1bmDIlbIcgkktU7iI1dMgh4Y1RTz8dth+++OIwT/+LX4S19CK5QOUuspsOPzzsTvnEE9CpE1xwAXToALfcAlu2xE4nhU7lLlJLRx0Fjz0Wdqk88EAYMSLcjL3jDvjss9jppFCp3EUywAz69Qsblf3lL7DffjB8OHTpAr/7XThhSqQuqdxFMsgsLJvcvgXxV74CZ5wRVtzMnBm2PRCpCyp3kSwwC1sML14Ms2eHA7+HDg17zM+erZKX7Eur3M1sgJmtNLPVZjZmF9cNNjM3s2o3tREpBPXqwf/8D7zwQhi5f/45DB4cVtyUloYzYUWyodpyN7P6wAxgINAVGGpmXau47ivABcDCTIcUyXf16oX9apYuhd//Prz5qaQkbF725z+r5CXz0hm59wFWu/sad98CzARKqrhuEnA9oJW+IjtRvz788Idh35o77giHiQwaBEccEVbbqOQlU9Ip9xbA+gqPN6Se+y8z6wW0cve5GcwmklgNGsBZZ4UdKG+5Bdavh2OPDdsbPPlk7HSSBOmUe1XHCv93fGFm9YCpwMXVfiGzc8yszMzKNm/enH5KkYRq1AjOPRdWrw7vcH3lFfjOd6B/f3jmmdjpJJ+lU+4bgFYVHrcE3qzw+CtAN+BxM1sLHAaUVnVT1d1vdfcidy9q3rz57qcWSZjGjWHkyHBoyJQp8OKL8K1vhSmbRYtip5N8lE65LwI6mFkbM2sEDAFKt3/S3d9392bu3trdWwPPAsXurjP0RGpojz1g9GhYsybsI/+Pf4SbriUlYVdKkXRVW+7uvhUYCcwHVgD3uvsyM5toZsXZDihSiL78ZbjsslDykyaFefheveDQQ+H66+G112InlFynA7JF8sB774XVNffeG0bzENbKn3JK+GjbNm4+qTs6IFskQZo0CVsLL1wYRu033BCWVY4ZE7YbPuSQcHjIq6/GTiq5QuUukmdat4ZLLglFv3YtTJ4ctjcYOxbat1fRS6ByF8ljBx0URvTPPlt10ffuDT//eVhqKYVF5S6SEBWL/vXX4cYbwxLL//3fcIhIr15wzTVhLb0kn26oiiTcunVhJ8r77it/Y1TPnuU3Yzt0iJtPaibdG6oqd5ECsn49zJq1Y9H36FFe9B07xs0n1VO5i8gurV9fPqJ/+unwXPfu5UXfqVPcfFI1lbuIpG3DhvKiX7AgPKeiz00qdxHZLVUV/Te+UV70nTvHzVfoVO4iUmtvvLFj0btDt27lRd+lS+yEhUfvUBWRWmvRAi64AJ56Kozob7oJmjaFCROga9dQ9JMmwUcfxU4qlancRSQtBxwAo0aFTcy2F/2++8L48eF0KR36nVtU7iJSYxWLfvp0eOABGDcudiqpqEHsACKS30aNgmXLwjYHXbrA6afHTiSgkbuI1JIZ3HxzOP91+PDyNfMSl8pdRGqtYcPwztdWreCkk8LeNhKXyl1EMuKrX4W5c+HTT6G4GP7979iJCpvKXUQypnPncFrUsmVw2mlaQROTyl1EMuq734WpU6G0NGw3LHFotYyIZNzIkbB8OVx3XVhBc8YZsRMVHo3cRSTjzMKbnPr1g3POKd+jRuqOyl1EsqJhw7AnzYEHhhU0a9fGTlRYVO4ikjX77htW0GzZElbQfPhh7ESFQ+UuIlnVqVMYwS9fDj/4AXz+eexEhUHlLiJZ179/2INm7lwYOzZ2msKg1TIiUifOPz+M3m+4IWwXfOaZsRMlm0buIlJnpk2DY44JK2j+/vfYaZJN5S4idWb7Cpo2bbSCJttU7iJSp5o2hTlzYOtWOOEE+OCD2ImSKa1yN7MBZrbSzFab2ZgqPn+emb1kZkvM7O9m1jXzUUUkKTp2DCP4FSu0giZbqi13M6sPzAAGAl2BoVWU993u/g137wlcD0zJeFIRSZRjj4Vf/ALmzYMxXxgySm2ls1qmD7Da3dcAmNlMoARYvv0Cd6/4g9WXAc9kSBFJphEjwg6SkyeHFTQ//nHsRMmRTrm3ANZXeLwB+Gbli8zsfOAioBHQLyPpRCTxpk2DlSvh3HOhfXs48sjYiZIhnTl3q+K5L4zM3X2Gu7cDLgeqPCrXzM4xszIzK9u8eXPNkopIIjVoEPaA376CZs2a2ImSIZ1y3wC0qvC4JfDmLq6fCZxY1Sfc/VZ3L3L3oubNm6efUkQSrWnT8O7Vbdu0giZT0in3RUAHM2tjZo2AIUBpxQvMrEOFh8cDr2QuoogUgg4dwjmsq1bB0KFaQVNb1Za7u28FRgLzgRXAve6+zMwmmllx6rKRZrbMzJYQ5t21Nb+I1Fi/fmEFzUMPwWWXxU6T39LaW8bdHwIeqvTc+Aq//2mGc4lIgTrvvLAHzZQpYQXNsGGxE+UnvUNVRHLOlCnhLNYRI+CJJ2KnyU8qdxHJOQ0awB//CG3bwsknawXN7lC5i0hOatJkxxU0778fO1F+UbmLSM5q3x5mz9YKmt2hcheRnNa3L9x8M/z5z3DppbHT5A+dxCQiOe/cc8MKmqlToUsXOPvs2Ilyn0buIpIXbrwRjjsOfvITePzx2Glqzh2WLIGrroKXXsr+62nkLiJ5YfsKmsMPDyto/vEPaNcudqpd27IFnnwSHnwQSkth3Towg/32g298I7uvrZG7iOSNffYJpzhB7q6gee89uOeecAO4eXPo3x/uuAN69gy/vvVWWL+fbRq5i0headcurKDp3x9OPTUsl2wQuclefz2MzB98MLzpauvWMDo/5RQoLg4Hk+y5Z91mUrmLSN45+mj41a/CjdVLLgl7wtcld3juufLplhdeCM936QIXXwwlJdCnD9SvX7e5KlK5i0heGj48nOI0bVrYg+acc7L7ep9+Gm7kbi/0N96AevXgiCPCSVLFxWFny1yhcheRvHXDDeEUp/PPD8Xat29mv/6774YdKktL4S9/gQ8/DNMrxx0XRufHHw/NmmX2NTNF5S4ieatBg3DzsuIKmvbta/c116wpnz9/6qnwrtivfz3cIC0uhmOOgS99KTP5s0nlLiJ5bfsKmm9+M6ygeeaZsC9NurZtg7Ky8umWpUvD8926weWXhxF6UVGYgsknKncRyXvt2sH994dVKaeeCvPm7XoFzSefwKOPhkKfMwc2bgw3P488MrwL9oQTcn8NfXVU7iKSCEcdFVbQDB8OF10EN9204+fffjuUfmkpzJ8PH30Ee+0FAweG6ZZBg2DffeNkzwaVu4gkxrBhO57idMwx5fPnCxaEKZgWLeBHPwrTLUcfDY0bx06dHSp3EUmU66+Hl1/e8V2gPXrAz34WCr1377AFQNKp3EUkUerXDytoxo6Fzp3D/Hnr1rFT1T2Vu4gkzt57w4wZsVPElWeLe0REJB0qdxGRBFK5i4gkkMpdRCSBVO4iIgmkchcRSSCVu4hIAqncRUQSyNw9zgubbQZe383/vBnwdgbjZIpy1Yxy1VyuZlOumqlNroPcvXl1F0Ur99owszJ3L4qdozLlqhnlqrlczaZcNVMXuTQtIyKSQCp3EZEEytdyvzV2gJ1QrppRrprL1WzKVTNZz5WXc+4iIrJr+TpyFxGRXcircjezVmb2mJmtMLNlZvbT2JkqMrP6Zva8mc2NnWU7M2tiZrPM7OXU/7fDY2cCMLPRqT/DpWZ2j5l9KVKO/2dmm8xsaYXn9jWzR8zsldSvTXMk1w2pP8cXzexPZtYkF3JV+NwlZuZm1ixXcpnZKDNbmfpeuz4XcplZTzN71syWmFmZmfXJxmvnVbkDW4GL3b0LcBhwvpl1jZypop8CK2KHqGQ68Bd37wz0IAfymVkL4AKgyN27AfWBIZHi/BYYUOm5McDf3L0D8LfU47r2W76Y6xGgm7t3B1YBY+s6FFXnwsxaAf2BdXUdKOW3VMplZn2BEqC7ux8MTM6FXMD1wFXu3hMYn3qccXlV7u6+0d2fS/3+Q0JRtYibKjCzlsDxwO2xs2xnZnsDRwF3ALj7Fnd/L26q/2oA7GFmDYA9gTdjhHD3J4F3Kz1dAtyZ+v2dwIl1Goqqc7n7w+6+NfXwWaBlLuRKmQpcBkS5ibeTXCOAa93909Q1m3IklwN7p36/D1n63s+rcq/IzFoDvYCFcZP81zTCN/e22EEqaAtsBn6Tmi663cy+HDuUu79BGEWtAzYC77v7w3FT7eBr7r4RwoAC2C9ynqqcBfw5dggAMysG3nD3F2JnqaQjcKSZLTSzJ8zs0NiBUi4EbjCz9YS/B1n5CSwvy93M9gJmAxe6+wc5kOd7wCZ3Xxw7SyUNgN7Ar9y9F/ARcaYYdpCawy4B2gAHAF82sx/GTZU/zOxnhCnKu3Igy57AzwjTC7mmAdCUMIV7KXCvmVncSED4iWK0u7cCRpP6yTrT8q7czawhodjvcvf7Y+dJOQIoNrO1wEygn5n9IW4kADYAG9x9+083swhlH9uxwGvuvtndPwPuB74VOVNF/zSz/QFSv9b5j/M7Y2ZnAN8DTvPcWMfcjvCP9Aup7/+WwHNm9vWoqYINwP0e/IPwU3Wd3+ytwhmE73mA+wDdUE39q3sHsMLdp8TOs527j3X3lu7emnBj8FF3jz4Sdfe3gPVm1in11DHA8oiRtlsHHGZme6b+TI8hB270VlBK+AtI6tcHI2b5LzMbAFwOFLv7x7HzALj7S+6+n7u3Tn3/bwB6p773YnsA6AdgZh2BRuTGJmJvAt9J/b4f8EpWXsXd8+YD+DbhZsSLwJLUx6DYuSplPBqYGztHhTw9gbLU/7MHgKaxM6VyXQW8DCwFfg80jpTjHsK8/2eEYhoGfJWwSuaV1K/75kiu1cD6Ct/7t+RCrkqfXws0y4VchDL/Q+p77DmgX47k+jawGHiBcM/wkGy8tt6hKiKSQHk1LSMiIulRuYuIJJDKXUQkgVTuIiIJpHIXEUkglbuISAKp3EVEEkjlLiKSQP8fVykw9QSpZMoAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(Ks2, np.array(scs), 'b-')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "显然，在K=2，即聚类数目为2时聚类效果最好。当然这也只是针对K-means的聚类情况。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "下面尝试先进行PCA降维，再聚类。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "K:2,CH_score: 3458.385024032641, Silhouette_Coefficient: 0.7039621761664018,time elaps:11\n",
      "K:4,CH_score: 2197.5133065970663, Silhouette_Coefficient: 0.5270905765051502,time elaps:10\n",
      "K:6,CH_score: 1782.6343217125957, Silhouette_Coefficient: 0.49620327401676995,time elaps:10\n",
      "K:2,CH_score: 3486.757100220883, Silhouette_Coefficient: 0.8623708745440023,time elaps:10\n",
      "K:4,CH_score: 2720.1424593027205, Silhouette_Coefficient: 0.5693743566113783,time elaps:10\n",
      "K:6,CH_score: 1685.6000336366465, Silhouette_Coefficient: 0.46105099543035427,time elaps:10\n",
      "K:2,CH_score: 3529.0356164862687, Silhouette_Coefficient: 0.7218452395565128,time elaps:10\n",
      "K:4,CH_score: 1732.050422699013, Silhouette_Coefficient: 0.42999083717887066,time elaps:10\n",
      "K:6,CH_score: 1839.1882090764946, Silhouette_Coefficient: 0.48719975182191166,time elaps:10\n",
      "K:2,CH_score: 3317.912118117484, Silhouette_Coefficient: 0.6727361367490524,time elaps:10\n",
      "K:4,CH_score: 2405.126679181873, Silhouette_Coefficient: 0.5470138767075347,time elaps:10\n",
      "K:6,CH_score: 1852.425096038561, Silhouette_Coefficient: 0.5088938838788698,time elaps:10\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD8CAYAAACMwORRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzs3XdclWX/wPHPfQYgG5TlwIWLSi1NbZgDB2oOKJ+yadlTqS3NxtP6tWxr61F7SptWpiVqauYuR0PNvXEgKCjK3mdcvz9uBBRQROAc4Pt+vXh5xn3u+8vx5nuuc93X9b00pRRCCCHqFoOjAxBCCFH1JLkLIUQdJMldCCHqIEnuQghRB0lyF0KIOkiSuxBC1EGS3IUQog6S5C6EEHWQJHchhKiDTI46cKNGjVSLFi0cdXghhKiVtmzZclopFXCx7RyW3Fu0aMHmzZsddXghhKiVNE2Lq8h20i0jhBB1kCR3IYSogyS5CyFEHeSwPvey2Gw2UlJSsFgsjg5FCCEcymw24+/vj9ForNTrnSq5p6Sk4ObmRqNGjbDZbNhsNkeHJIQQNU4pRU5ODsnJyQQHB1dqH06V3C0WCw0bNiQpKYm4OP2CsKZpDo5KCCFqnlKKrKwsTp06RceOHS/59U6V3AGys7M5evQo7u7uGAxySUAIUX8VFBSwevVqQkJCCAi46ND2czhd9iwoKEDTtPITuyUXUo+CLA8ohKjjDAYDmqaRmZl56a+thnguy0XXdM05A+kJkLwX7NInL4QQZXG6bpmL8mkKmgFSDsPJXRAYDkazo6OqEzZt2oTZbKZz586V3kd2djajR48uun/y5EmGDBnCM888Q0FBAc8//zx79uzBx8eHd999lyZNmlRB5EJUTlWc8wDLli3js88+w26307NnTyZOnAjg0HPe6VruFeLdGALaQ0E2JG7Xu2rEZdu0aRPbtm27rH14eHgwb968op+QkBAiIiIAmD9/Pt7e3ixZsoS7776bDz74oCrCFqLSquKcT0tLY+rUqXz22WfExMRw5swZ/vzzT8Cx57zTttw/25zK4dSCC29kbw6WHCABzB5guPB40FZ+Lvy7q98Ftzl+/Djjxo3j6quvZvv27QQGBvLhhx/i5uZWattjx47x2muvkZqaisFgYMqUKTRt2pSpU6eyfv16NE3jwQcfJDIykk2bNjF9+nQaNmzI/v37iYiIoE2bNnz77bfk5eXx4Ycf0qxZM1544QVcXV2JjY0lJSWFSZMm0atXL/Lz83n99dfZvXs3JpOJSZMm0a1bNxYuXMjatWvJzc0lISGBvn37FrUaNm7cyPTp0ykoKKBZs2a89tpruLu7ExkZydChQ/ntt9+wWq1MmTIFFxcX5s2bh9FoZMmSJTz77LOcOXOGGTNmYDQa8fT05Msvv7zw/8d54uLiSElJoUuXLgCsXbuWsWPHAtC/f3/efPNNlFIyIqrQ1/FfE5dbobIhFda8QXPuaXbPBbeRc/7yzvmEhASaN2+Ov78/AD169GDlypX06NHDoee80yb3CjGYwMVTb8EXZINLAzBcfhfNsWPHePvtt3n55ZeZNGkSK1eu5Oabby613bPPPsuYMWOIiIggPz8fu93OypUr2b9/Pz/++CNpaWmMGjWqKLkdOHCABQsW4OPjw+DBg4mKiuK7775j9uzZfPfddzzzzDOA/sf2xRdfEB8fz5gxY+jRowdz5swB9JbAkSNHeOihh/j5558B2LdvH3PnzsXFxYVhw4Zxxx134Orqyqeffsqnn36Ku7s7n3/+OV9//TUPP/wwAH5+fsydO5c5c+bw5Zdf8sorrzBy5Ejc3d2LulWio6P55JNPCAoKIiMjA4AjR47w9NNPl/m+zZo1C29v76L7v/zyCwMHDiw6kU+ePElQUBAAJpMJT09P0tLS8PO78AeuqH5yzo8GKnfOh4aGcuTIEY4fP05QUBCrV68umojpyHPeaZP7xVrY57BZ4ORuKDgJ/q3BO+Syjt2kSRPat28PQHh4OCdOnCi1TXZ2NqdOnSrqcnB1dQVg69atDBo0CKPRSMOGDenatSu7du3C09OTK664omg4U9OmTbn++usBaNOmDZs2bSra98CBAzEYDDRv3pymTZty5MgRtm7dyqhRowBo2bIlISEhRXMBunfvjpeXFwCtWrXixIkTZGZmcvjwYe69915An0PQqVOnomOcjTs8PJxVq1aV+T5cffXVvPjiiwwYMIB+/foVHXvevHkVeh+XLVvGG2+8ccFtpNVe7GIt7Ook57yusuf8Cy+8wFNPPYXBYKBz584kJCSUu21NnfNOm9wvidEMwVdB8n5IOQS2PPBtAZV8E83m4ta/wWDAarWW2qa8UT0XGu3j4uJyzn7P3jcYDOfMxj3/P1/TtArv12g0YrPZUErRo0cP3nnnnQu+5uz2ZXnxxRfZsWMH69atY+TIkcybN4/U1NQKtdz379+PzWYjPDy86PmgoCBOnjxJcHAwVquVrKwsfHx8yv29RM2Rc15X2XO+d+/e9O7dG4Aff/yxaCi3I8/52nlBtSwGIwR2AK9gSD8Op/eDslfb4Tw9PYu+goF+VTw3N5cuXbrw66+/FtXJ2bJlC1ddddUl7Xv58uXY7Xbi4+NJSEigRYsWdOnShaVLlwJw9OhRkpKSuNBiJx07dmTbtm0cO3YMgNzcXI4ePXrB43p4eJCdnV10Pz4+no4dOzJ+/Hj8/PxISkoqasWU9XN+l0xkZOQ5++/duzeLFi0CYMWKFXTr1k1a7rWInPPln/NnzpwBICMjgx9++IHo6GjAsed83Wi5n6Vp0DAMTK6QGqd31wR20Pvmq8Ebb7zBq6++yrRp0zCZTEyZMoWIiAi2b9/OrbfeiqZpTJgwgUaNGnHkyJEK77dFixbcd999pKSk8OKLL+Lq6sptt93Ga6+9RnR0NCaTiddee+2c1sv5/P39ee2114qGIAI88sgjF/zj6NWrF08++SRr167l2WefZfbs2cTFxaGUonv37rRr167Cv8Ovv/7K9OnTz3ksKiqK5557jiFDhuDj41NuC0s4Lznny/b2229z4MABAB566KGiYzrynNcuOmmomnTt2lWdvxLTiRMncHFxITY2Fg8Pj8s7QNYpOH0QzA0g6Ao94dcCL7zwAjfddBMDBgxwdChC1Ag558uXnp7Oli1bGDp0KK1atQJA07QtSqmuF3tt3emWOZ9noJ7Ubfn6WPiCLEdHJIQQNaZudcucr4EvBHeEk3sgcScEtocGlz4EafLkyaUmOtx5552MGDGiqiIt8vrrr1f5PoW4VHLO1351O7kDuHhASEc4tUdP8o3CwDPoknbx/PPPV1NwQjgnOedrv7rbLVOSyVUfKunmo/fDpx1zdERCCFGt6kdyB33ETFC43hefdkxP8lI2WAhRR9X9bpmSNAM0agtGV0iP1y+2BnS4aE0aIYSobepPy70kv+b6ePjcdEjaAbaLFCgTQohapn4md9BnsgZ2AEteYdngHEdH5HBVUf4UYOnSpURHR3PLLbfw8MMPk5qaCuhjdh988EFuvvlmHnzwwaLCTEI4SlWd8x999BH9+/ene/fu5zxeUFDAU089xZAhQ7jjjjs4fvx40XMzZ85kyJAhDB06lA0bNlx2DOerv8kdwN1fv9Cq7JC4A/LSHR2RQ1XFiW61Wnn77beZNWsWP/30E23btuX7778H9Doc3bt3Z/HixXTv3p1Zs2ZVRdhCVFpVJfdevXrx3XfflXq8vHruhw4dYtmyZcTExDBjxgwmT55cbr2bynLaPveQ7R/ilnawSveZ59uGxE6Pn/ugqyeEdNKrSp7cxXGLN+OefEFqW1eytvXZGc+5ubn4+vqSnZ1NaGgoAGvWrOHzzz8HYNiwYdx///1MmDDhMv5H6xbDF1+iXaQWyqVSLVpgv2/0BbeReu6Xv4ZByeqTJZVXz33NmjVERkbi4uJC06ZNCQ0NZdeuXeXupzKcNrnXKJObPtnp1F5IOsyxY3FS25rK13N//vnnueWWW2jQoAGhoaE899xzAKSkpBSVfw0ICCAlJaUq/vdEFZB67qOBy1/D4Hzl1XM/deoUHTt2LNrubPXIquS0yb1UC7u6Gc0QfCUkp9IksCHtA8yglNS2vsTa1haLhblz5zJ37lyaNm3Km2++yaxZs3jwwQcv+PYLLtrCrk5Sz113uWsYVER55Yyrulqk0yZ3h9AM4N8as2sDyEwEWwEGgya1rS+htvXZcqvNmjUDYMCAAUVdMf7+/iQnJxMQEEBycnLRsmTC8aSeu+5y1jAoS3n13IOCgkhKSira7uTJkwQGBpa7n8qo3xdUy6JpeivevxXknIGMxDLrwktt67JrWwcGBnL48OGiLpc///yzqJpdydrWixYtok+fPpf0vgjHknP+wmsYlKW8eu69e/dm2bJlFBQUkJCQQFxcHFdeeWUF3qmKk5Z7ebwbg9EFbCsh5yRYcvXywSVIbevSAgMDefjhh7nvvvswmUyEhIQUFYYaM2YMkyZNIiYmhuDgYKZMmVLh90Q4BznnyzZ16lSWLl1KXl4e/fr1Izo6mnHjxpVbzz0sLIwBAwYwYsQIjEYjzz33HEZj1U6mrFA9d03TIoEPASMwUyn11nnPhwJfAb6F2zyrlFp6oX1Wez33qpKXoRcdQ9PLF7h6VevhpLa1qG/knC9ftdZz1zTNCEwDBgHhwChN08LP2+wFYK5S6mrgdmA6dYWbtz5U0mCEpJ16V40QQji5inTLdANilVKHATRNmwMMB/aU2EYBZzuffIDSl9prscnvTGXb1q1gzdP7340u3HnPaKltLeosqede+1UkuTcB4kvcTwC6n7fNy8ByTdMeBTyAfpUNyFHL/l1IUW1ruw2S90FuKvg01atKygLPog6Seu6Op5S6rHxYkdEyZWWv8484CvhSKdUUGAx8o2laqX1rmvagpmmbNU3bnJycXGqnZrOZvLy8CoTkIAYjBIbrdWnSE+D0gTJH0gghxOVQSmG328nNza30PirSck8AmpW435TS3S5jgMjCoP7QNM0NaAScOi/gT4FPQb+gev6B/P39OXjwIJmZmWWOs3UapgBwAVKSICsX/EJBk7LBQoiqoZQiNze3aEixyXTpAxsr8opNQBtN01oCx9EvmN5x3jbHgAjgS03TOgBuQOmm+UUYjUaaNWvGhg0byMzMLJoB57TOxMLh36GBD7QZqNepEUKIKpKfn4+fn1+lJjhdNLkrpayapj0C/Io+zPFzpdRuTdNeBTYrpRYBTwKfaZo2Ab3LZrSqZGeRh4cHI0eOZOfOnWRlZVVmFzWnVSto0QI2z4L4edD9Ib0vXgghqoCXlxcdO3Yss4jbxVRonHt1KGuce611cjfMvhXyM+G2r6F1X0dHJISoo6psnLuogKAr4IGV4BsK346EbaXrOgshRE2S5F5VfJrA/b9A8xtgwVj47R1ZgFsI4TCS3KuSmw/c+SN0vB3WTIafHwObxdFRCSHqISkcVtVMLhD1iX5hdd17kHECRn4lI2mEEDVKWu7VQdMg4kW4+QM4tAa+HAyZVbvKihBCXIgk9+rU9T4Y9T2cPggz+0HyfkdHJISoJyS5V7e2A2H0ErDmwqwBELfR0REJIeoBSe41ock1MGYFeDSCr4fD7hhHRySEqOMkudcU/5Z6gm98DcwbDRv/K0MlhRDVRpJ7TXL3h3sWQIdhsPx5WPasXkZYCCGqmCT3mmZuoA+N7DEe/voE5t6jr88qhBBVSJK7IxgMEPkGDHwT9i2Br4ZBtizfJ4SoOpLcHem6cTDyS0jcDrP6Q8phR0ckhKgjJLk72hUj4J6FkJsCM/vD8S2OjkgIUQdIcncGza/TR9K4uMOXN8P+ZY6OSAhRy0lydxaN2sCYldCoLcwZBZtmOToiIUQtJsndmXgF6bNZw/rBkomw8mWwywLcQohLJ8nd2bh6wu3fQ5fRsP59iHkIrAWOjkoIUctIyV9nZDTpFSV9msHq1yAzEW6bDQ18HR2ZEKKWkJa7s9I0uGkSRP0Pjv0BXwyC9ARHRyWEqCUkuTu7TrfrqzulxetDJZN2OToiIUQtIMm9NmjdB+4vHB75eaS+AIgQQlyAJPfaIvhKeGAl+DaDb2+Fbd87OiIhhBOT5F6b+DTRW/DNr4cFD8Nv70rZYCFEmSS51zZuPnDnT9DxNljzOvz8GNisjo5KCOFkZChkbWRy0UfR+DSFdVMgI1EvQObq6ejIhBBOQlrutZWmQcRLcPP7cGgVfDkYMk86OiohhJOQ5F7bdb1fn9F6+iDM6gfJBxwdkRDCCUhyrwvaRcLoxfqKTrP6Q9wfjo5ICOFgtS6551hyWJewztFhOJ8mXfSywR6N4OvhsDvG0REJIRyo1iX3mTtnMm7VOL7c9SVKhgGey7+lnuAbd4Z598Ef0xwdkRDCQWpdcn+o00MMaD6AKVum8Obfb2Kz2xwdknNx99dXduowFH59Dn55FuQ9EqLeqXXJ3dXoyru93uXe8Hv5ft/3TFg7gVxrrqPDci7mBvrQyO5j4a8ZMG+03h8vhKg3al1yBzBoBiZdO4lnuz3L2vi1jPl1DGdyzzg6LOdiMMKgt2DgG7B3kd4Pn5Pi6KiEEDWkVib3s+7scCfv93mfA6kHuGvpXRxNP+rokJzPdeP1VvyJbfpImpQjjo5ICFEDanVyB4gIjWDWwFlkW7K5+5e72XZqm6NDcj5XROn98Nmn9QR/fIujIxJCVLNan9wBOgV0Yvbg2Xi7ePPA8gdYEbfC0SE5n+bX6SNpzA3gy5th/zJHRySEqEYVSu6apkVqmrZf07RYTdOeLeP59zVN21b4c0DTtLSqD1WXU2DlRFrpi4Oh3qHMHjyb9v7teXLtk3yz55vqCqH2CmgLY1ZCo7YwZxRs/tzREQkhqslFk7umaUZgGjAICAdGaZoWXnIbpdQEpVRnpVRn4GNgfnUEC/DFhqP0fm8tk5fsITX73IWj/dz8mDlgJn1D+/LOpnd4+++3Zajk+byCYPQSaB0BiyfAylekbLAQdVBFWu7dgFil1GGlVAEwBxh+ge1HAdW2ksTwzo0Z1qkxs9Yf4aZ31vDxqoNk5xeXvHUzuTGl1xTu6nAXs/fOZtJvk8iz5lVXOLWTqyeMmgPX3APrp0LMQ2AtuPjrhBC1RkWSexMgvsT9hMLHStE0rTnQElhdzvMPapq2WdO0zcnJyZcaKwBN/dx5b2Qnlj1xE9e1bsiUFQfo9e4avtp4lAKrHQCjwcgz3Z7h6WufZtWxVTyw/AFS81Irdbw6y2iCoR9Bnxdgxw/66k556Y6OSghRRSqS3LUyHivve/ztwI9KqTL7QpRSnyqluiqlugYEBFQ0xjK1DfLi03u6Mn/c9bQO8OT/Fu0mYupaYrYmYLPr4d0dfjdTek9hX8o+7lp6F8cyjl3WMescTYNeT8GITyBuA3w+CNKPOzoqIUQVqEhyTwCalbjfFDhRzra3U41dMmW5JtSPOQ/24Kv7u+HtZmbCD9sZ8tE6Vu09iVKK/s37M3PATDIKMrhr6V3sSN5Rk+HVDp1HwZ3zIO0YzOwHSbscHZEQ4jJVJLlvAtpomtZS0zQX9AS+6PyNNE1rB/gBNV5vVtM0erUN4OdHbuTjUVeTZ7Ex5qvNjPzkDzYdTaFzYGe+GfQNHmYPxvw6hlXHVtV0iM6vdV+4/xdAwReD4PBaR0ckhLgMF03uSikr8AjwK7AXmKuU2q1p2quapg0rsekoYI5yYKlGg0FjaKfGrJjYi8lRV3IsJYeRn/zB/V9uIjfHn9mDZ9PGrw0T1kzg273fOipM5xV8FTywUl++b/YtsH2OoyMSQlSS5qhc3LVrV7V58+ZqPUZugY0vNh7hk7WHyMy3MrxTY8b1bc5/d73Mmvg13Bt+LxO7TsSg1Ym5XFUnNw1+uAuOroO+L0LPJ/X+eSGEw2matkUp1fWi29Xl5H5Weo6FT34/xBcbjmCzK26/timq4QIWHp7HgOYDeKPnG7gaXWskllrDmg8LH4Gdc6HLaBg8RR9hI4RwqIom93rx1+rjbuaZyPaMvr4FH606yHd/x+NivJYeVzdgedzXJOcm81Gfj/B183V0qM7D5ApR/9O7aNZPhYxEuPVzfYy8EMLp1av+iCBvNyZHXcXKib3oFx7Mmr/D0ZLvZvupndz1y93EZ8ZffCf1icEA/f4PhkyF2BXw5RDIOuXoqIQQFVCvkvtZLRt58PGoq1n86I109u9N1tExxKWdYuTCO9h+UoZKlnLtGLj9Ozh9QB8qefqgoyMSQlxEvUzuZ13ZxIev7+/GN3fdTrPcp8nM1bjrl/uYsm6+rM96vnaDYPRiKMjWywYf+9PREQkhLqDWJffqSLrXt27EkrHRvHLt/zDZgvni0Mv0+fQNNsServJj1WpNusADK6CBP3w1DPYsdHREQohy1LrknrF0KUf+dRup33+PLb3qaqFomsatnTuw7q65tPW6ljNuc7hv4cvcOfMPtsdXWwXj2se/lV4XPqQTzL0X/pju6IiEEGWodcnd4OqKyssj6ZVXOdjzJo5PfJKsdetRtqop7evp6sHcqP9xS9hIXBv9xs6CGQyf/hvjvt3CoeSsKjlGrefREO5dBO2HwK//gWX/Abvd0VEJIUqolePclVLk7dlD+vwYMhYvxpaejik4GJ/hw/GNGoFLixaXHZ9Sis93fc4H/3xAsEs4J/bfRl6BGyO7NOXxfm0I8Wlw2ceo9ew2+PU5+OsTCB8OUZ+C2c3RUQlRp9WbSUz2ggKyVq8hLWY+2evWg91Og2uuwTc6Cq/IQRg9PS5r/0sPL+WFDS/Q2KMp4YaJLNiUCxqMvr4FY3u1xs/D5bJ/h1pNKfhjGix/Hpr1gFHfg7u/o6MSos6qN8m9JMupU2QsWkTa/BgKDh9Ga9AA7wH98YmKxr3btWiGyvVCbUraxONrHsfV6MpLXaeweLOR+VsT8HQx8VCvVtx3Q0s8XOvFfLDy7ZqvL/rh2xzu+hH8Wjg6IiHqpHqZ3M9SSpG3Ywdp82PIWLIEe1YW5iZN8BkxAp+oEbg0bXrJ+zyUdoixK8eSlp/GlF5TCDR15r3l+1mx5ySNPF15LCKM268NxcVU6y5jVJ24jfD9KDCa4Y650OQaR0ckRJ1Tr5N7Sfa8PDJXrCQ9JobsP/4ApXDv3h2fqBF4DxiAwd29wvtKzklm/KrxHEg9wAs9XuDWtreyJS6Vt5ft4+8jKTTzb8CT/dsxrFNjDIZ6WmgreT/MvhVyTsPIr6DtAEdHJESdIsm9DJYTJ0hfuJC0mAVYjh3D4OGB16BIfKOiaHDNNWgVqHyYY8nhyd+eZP3x9fz7qn/z6NWPAvDbgWTeWbafPYkZtA/24unIdvRpF1ihfdY5mUnw3b/0RT9unqoXHhNCVAlJ7heglCJ3yxa922bZMlRODi7Nm+MTFYXPiOGYg4Mv+Hqr3crrf77OTwd/4uZWN/Pq9a9iNpqx2xWLdyYyZfl+4s7kcG0LP56ObM+1LerhBcb8LJh3L8SuhJ6ToO8LUjZYiCogyb2C7NnZZCxfQfr8+eRs2gSahsf11+MTHYVXv34YXMsuBayUYubOmXy09SO6BXfj/T7v4+3iDYDFZueHTfF8uOogyZn5RLQPZNLAdnQI8a7JX83xbBZYPAG2fgMdb4dhH4Opno8uEuIySXKvhIL4eNJjFpC2IAbriUQM3t54Dx6Eb3Q0blddVWYXy8+HfualjS/RwrsFM/rNINijuNWfU2Dly41HmbH2EFn5VkZ0bsLE/m1p5l/xfv5aTyn4/V1YMxla9oLbvgE3H0dHJUStJcn9Mii7nZy//iItJobM5StQeXm4hLXGNyoKn2HDMAUEnLP9X4l/8cSaJ3A3uTOt3zTa+7c/5/m0nAI++e0wX2w4gl0p7ugWyiN92xDgVY8WCNn2HSx6FBq10xfj9mni6IiEqJUkuVcRW2YmGb/8QnrMAnK3bgWjEc+ePfGJisKrT280F72b4WDqQcauHEtmQSZTe0/lhiY3lNpXUnoeH60+yA+b4nE1GRhzY0v+fVMrvN3MNf1rOcah1fDDPeDqpY+FD7rC0REJUetIcq8G+YePkB4TQ/rChVhPncLo64v30KH4Rkfh1qEDJ7NPMn7VeGLTYvm/6/6PqDZRZe7ncHIWU1ccYPGORHzdzYzvHcbd1zXHzWys4d/IARJ3wLcjwZIDt82GVr0cHZEQtYok92qkbDayN24kbf58slauQlksuLZvj290FMaBfXl6x2tsPLGRhzs9zLhO48odDrnreDrv/Lqf3w8kE+LjxhP92nDLNU0xGev4RKi0eD3Bn4mF4dOg022OjkiIWkOSew2xpaWRvnQp6fNjyNu1C8xmPHrdxNIOuXzS4C+GtB3Oy9e9jNlYftfLxkOneWfZfrbFp9EqwIOnBrQj8srguj1GPjcNfrgLjq6DiJfgxokyVFKICpDk7gB5Bw6QHrOA9EWLsJ05Q4GPO7+2yyUlojMvjPoUT5fyF5dWSrF8z0ne/XU/saey6NTUh2ci23N9WKMa/A1qmDUfFo6HnfOgy30w+D0w1vMaPUJchCR3B1IWC1nr1pEeE0PG6tVoNjsJzdxoe+fDNI26HaNP+UMBbXbF/H8SeH/FAU6k59GzTSOeGtiOjk19a/A3qEF2O6x+Fda/D20j4dbPweXyKnkKUZdJcncS1pQUts/+mNM/zSX0pB1czHj364dPVDQe11+HZiz7ImqexcbsP+OYtiaW1BwLQ64KYeKAtrQOKL/1X6ttmglLn9JXeLpjLngGOjoiIZySJHcnsz9lP29882+6bkmn7z4zhoxsTEFB+Awfjk/UCFxbtizzdZl5Fj5bd4SZ6w6Tb7Xzr65NeSyiji4Wsm8p/Hi/ntjvmg+NwhwdkRBOR5K7E0rKTmLcqnHEnz7MW8Zbabsxgax16/QFRq6+Gp/oKLwHDcLoWbp1fjorn/+ujuXbv+IwaJq+WEjv1vi617Hp/Amb4bvbQNlh1BwI7e7oiIRwKpLcnVRmQSYT1k7gr8S/GN95PGOCo8lY/LO+wMihQ2hubngN6I9vdDTu3bqVWmAkPiWHD1Ye1BcLcTXxcK/W3HdDC9xd6tCFyJTDMPsWyDgB0Z9B+DAT8QlDAAAgAElEQVRHRySE05Dk7sQsNgsv//Eyiw4tIrpNNC/0eAGTZtIXGImJIWPJUuyZmRdcYGR/Uibv/rqflXv1xUIejwjjtrq0WEj2afj+dr0lH/km9Bjr6IiEcAqS3J2cUorp26fzyfZPuKHxDUzpPQUPsz5KxJ6XR+bKVfoCIxs36guMdOumd9uct8DIlrgU3v5lP38fTSHU350nB7RlaMc6slhIQQ7M/zfsWww9xsOA16GSSyUKUVdIcq8l5h+cz6t/vEobvzZMi5hGoPu5o0QsiYmFC4zEYIk7hsHdXV9gJDq6aIERpRRrCxcL2Vu4WMgzke3p3S6g9k+Esttg2X/g7/9B+AiI+h+Y3RwdlRAOI8m9FtlwfAMT107E29WbGREzCPMrPUpEKUXuP/+QNn8+mb8sw56Tg7l5qF6pcvhwzCEh2O2Kn3ecYMryAxxLyaFbC3+ejmxH19q+WIhS8Md/YfkLEHod3P4duNfy30mISpLkXsvsPbOX8avGk2fN44M+H9AtpFu529pzcshYvpz0+THk/P23vsDIddfhEx2NV78IrCYXftgcz0eFi4X066AvFtI+uJYvFrLrJ4h5GPxawJ0/gl9zR0ckRI2T5F4LJWYlMm7VOI5mHOW1G17j5lY3X/Q1BfHxpC9YSHpMDJYTJzB4eeE9ZDC+UVHY24fz5cY4PvlNXywkqnMTJtT2xUKOboA5o8DoCnfOhcZXOzoiIWqUJPdaKqMggyfWPMGmpE08dvVjPHDVAxXqN1d2Ozl//62XPPh1ub7ASOvW+EaNgP6D+HRPBl9uOIpdKe7s3pzxfcJq72Ihp/bBt7dCTgqM/BLaDnB0RELUGEnutViBrYCXNr7EksNLuLXtrTzf/XlMhoqPY7dlZRUvMPLPP/oCIzfeiBo4hP9ZmzBnWxKuJgMPFC4W4lUbFwvJTNLLBp/cDTdPhS6jHR2REDVCknstp5Ti460f89nOz+jZpCfv9XoPd/Old6fkHzmid9ssWID15Em9aFm/SL7zv4ovT7ni5+HC+D5h3NWjFi4Wkp8J80ZD7Eq46Sno87yUDRZ1XpUmd03TIoEPASMwUyn1Vhnb/At4GVDAdqXUHRfapyT3ipl3YB6T/5xMW7+2TO83nUYNKlcCWF9g5A/SY+aTuXIVqqAA1SqMlc2vZZZbOzwCG/FEv7ZEX9Okdi0WYrPA4gmw9RvoNAqGfgSmOlaSQYgSqiy5a5pmBA4A/YEEYBMwSim1p8Q2bYC5QF+lVKqmaYFKqVMX2q8k94r7PeF3Jv02CT9XP2b0m0Er31aXtT9bejoZS5eSNj+GvJ07UUYje5p3ZF5AZ1Ku7MKTg8IZeEUtWixEKfjtHVj7BrTqDf/6Btxq+cggIcpRlcn9OuBlpdTAwvv/AVBKvVlim3eAA0qpmRUNUJL7pdl9ZjfjV46nwF7AR30+omvwRf9vKyT/4EHSYhaQvmghttNnyGjgxYrGV3OsewSj7+hbuxYL2Tobfn4cAtrDnfPAu7GjIxKiylVlcr8ViFRKPVB4/26gu1LqkRLbLEBv3d+A3nXzslJqWRn7ehB4ECA0NLRLXFxcxX8jQUJmAuNWjSMhM4HJN05mUMtBVbZvZbGQtX49qT/NJ3PNWgw2Kwd8m3KkS28ixt/NVeGhVXasahW7CubeA24++lj4oHBHRyRElarK5D4SGHhecu+mlHq0xDaLAQvwL6ApsA64UimVVt5+peVeOen56Ty2+jH+OfUPE7pM4L4r7qvy7hNraipnFi7i2Ldz8Yw/jMVg5Gj7rnS4bxRhg/uVu8CI00jcoY+kseTC7bOh5U2OjkiIKlPR5F6RK2cJQLMS95sCJ8rYZqFSyqKUOgLsB9pUNFhRcT6uPnw64FMiW0Ty/pb3mfzXZKx2a5Uew+TnR9Doe7l2xRIC58zleM9BBMXuwvrUE/xzXU+OvPEO+YePVOkxq1RIR3hgJXiHwDfRsGOuoyMSosZVpOVuQu9yiQCOo19QvUMptbvENpHoF1nv1TStEbAV6KyUOlPefqXlfnnsys4H/3zAF7u+oHfT3rx909uVGipZUadSslg44wdcViyly8l9GJUdc8dONLw1Gu/Bg8tcYMThclNhzl0Qtx4iXoIbJ8pQSVHrVfVQyMHAB+j96Z8rpSZrmvYqsFkptUjT+wWmAJGADZislJpzoX1Kcq8a3+/7nrf+fotw/3A+jvi40kMlKyo+JYdP5v9N3rIlDDy2mWYZSfoCI/374xsdhXv37qUWGHEoaz4sGKvXpel6Pwx6F4x1aGETUe/IJKZ6ZM2xNTz9+9M0bNCQGf1m0NKn7PVYq9K+pAzeW7afY39sZljiP/RK2IYpJwtT4xB8R4zAJyoKl2bNLr6jmmC3w6qXYcOH0HYQ3DoLXDwcHZUQlSLJvZ7ZmbyTR1Y/gk3Z+KjPR1wTdE2NHPfsYiFbD51kaNZB7kjZgceOLfoCI9dei090NN4D+mPwcIJk+vdn8MvTENIZ7pgLngGOjkiISybJvR6Kz4xn3MpxnMg6wRs932Bgi4E1clylFGv3J/P2sn3sS8qkh5eNieog/utWYImL0xcYiYzENzqKBl26OHZy1L4l8OMY8AqCO3+CRqVr5wvhzCS511NpeWk8uvpRtiVvY1LXSdwTfk+NJdPSi4X48Z/QAoI2riRz6S/6AiOhofhGjdAXGGnsoElGCZvhu3/pM1tHzYHQ7o6JQ4hKkORej+VZ83hu/XOsiFvBHe3v4Olrn8ZoqLmx6QVWOz9sOsaHq2I5nZVPvw5BTLoplJAdf5A2P4acv/4qXGCkBz5R0Xj174fBrYaXzjtzSC8bnHECbpkJHYbW7PGFqCRJ7vWcXdmZsnkKX+/5mr7N+vLWTW/RwNSgRmPIKbDyxYajfLL2EFkFVqKubsKEfm0JykklfcECfYGR48cxeHriPXgwvtFRuHXqVHPdNtmn4fvb9Zb8oLeh+0M1c1whLoMkdwHAt3u/5e2/3+aqRlfxccTH+LvV/NqjqdkFfPLbIb7cWLxYyCN9w2jobiZn02bS588nY/lyVG4uLq1a4RM1Ap9hwzEHBV5855erIAd+egD2L4HrHoH+r4EzDeUU4jyS3EWRVXGreGbdMwS6BzKj3wyaeztm7dHE9Fw+WnWQuZsT9MVCerbi3z1b4uVmxpaVTeayX0iLWUDuli1gMODR80Z8o6Lx7NsHg0s1lvG12+CXZ2DTZ3BFFIz4BMw13E0kRAVJchfn2J68nUdXPYpC8XHfj+kc2NlhsRxKzmLq8gMs2ZmIn7u51GIhBUeP6pUqFy7EmpSE0ccH75tvxic6Crfw8OrptlEKNn4EK16CZj2g97PQoqdMeBJOR5K7KOVYxjHGrhzLyZyTvNXzLfo17+fQeHYkpPHur/tZd/A0jX3ceKJ/W6KvLl4sRNlsZP/xJ+nz55O5ciWqoADXtm3xiY7CZ+hQTA0bVn1QO3/UywYXZEEDP2g/BDoMh1a9wFRL15wVdYokd1GmlLwUHl39KDuTd/L0tU9zV/hdjg6JDbGneWfZPrYnpBMW6MmkAe0YeEXQOS10W0ZG8QIjO3aAyYRnr174RkfhedNNaOYqXAe2IAcOrYI9i+DAMsjPAFdvaBsJ4cMhLALMNXtxWoizJLmLcuVac/nPuv+w6tgq7upwF09d+xQGzbEXEZVS/Lo7iXd/3c+h5Gw6NfPlmch2XN+6dK2c/NhY0mJiSF+0CFvyaYz+/vgMHYpPdDRu7dpWbWDWfDi8Vk/0+5foxcjMHtCmv57o2wwAVycsmibqLEnu4oJsdhvvbn6Xb/d+S//m/XnjxjdwMzn+IqLVZmf+P8d5f+UBEtPz6NmmEU8PbM9VTX1KbausVrLWrSM9ZgGZa9aAxYLbFVfo3TZDhmD09a3a4GwWOLpOT/T7FkN2MpjcoHWEnujbReqLhAhRjSS5iwr5evfXvLf5PToFdOKjvh/h5+bn6JAAyLPYmP1nHP9dE0tajoUhHUN4sn9bWgWU3Uq2pqaS8fNi0hbEkL9nL5rZjGdEBL7RUXhcfz2aqYovjNptcOwPPdHv/RkyT4DBrK/hGj5c76t3r/lhp6Luk+QuKmz50eX8Z91/CPEMYUbEDJp5O0k1RyAjz8LM3w8zc/0R8q12bru2GY9HtCHIu/xvGXn79pE2fz4ZPy/GlpqKKTAQn+HD8ImKwrXV5S0uXia7HY5vgT0LYO8iSDsGmhFa3FiY6G/Wa9kIUQUkuYtLsvXUVh5d/ShGzch/+/6XqwKucnRI50jOzGfamli+/SsOg6Zx3w0tGdurNT7u5V9IVQUFZP72G+nzY8j6/Xew2WjQqZNeqXLwIIxeXlUfqFKQuB32LNQT/ZlYQIPQ6/RE32Eo+DSp+uOKekOSu7hkR9OPMnblWE7nnuadm96hT2gfR4dUSnxKDu+vOEDMtuN4upp4uFdr7ruhBe4uF+52sZ4+Tfqin0mPmU/+wVg0V9fiBUZ69KieBUaUglN7ixP9qT364026Qvgw6DAM/Ku/9r6oWyS5i0o5k3uGR1Y9wp6UPTzb7VlGtR/l6JDKtC8pg/d+3c/KvacI8HLl8Yg23HZtM8zGCydppRR5u3aTHjOf9MVLsGdkFC8wMmIELqGh1Rf06YPFiT5xu/5YcEc90YePgEay7LC4OEnuotJyLDk8s+4Z1sav5b4r7uOJLk84fKhkeTYfTeHtZfvYdDSV5g3deXJAO26+KgSD4eKzWO35+WStWkVazAKyN2wAux33rl3xHj4M17AwzCEhmAIC0IzVUFEz9ah+IXbPQkjYpD8W0KEw0Q+HwHBZ71WUSZK7uCw2u423/n6LOfvnENkiktdvfB1Xo3PO0Dx/sZDwEG+ejmxHr7YBFS5VYDl5kvQFC0mPiaHg6NHiJ4xGTEGBmINDMIeEYA4JxhQcgrlxCObgYEwhIRh9fS+vJEL6cX1o5Z6FELcRUODfurjrpvHVkuhFEUnu4rIppfhq91dM2TKFawKv4aO+H+Hj6rzjuM9fLKR7S3+ejmxPl+YVH96plKLg8GEsx49jOZGIJSkRa2ISlsRELElJWBMTURbLOa/R3NwKE31w0YdA0e3CD4EKLzOYdao40R9ZB8oGPqHFLfomXaVqZT0nyV1UmWVHlvHc+udo4tmEGf1m0NSrqaNDuqDzFwvpHx7EUwPb0Tbo8kfHKLsdW0oKlsSkwsSfWHz7ROEHQHKyPjyyBIO3t97yP/shENK48FtAMObGjTEHBqKdX/kyJwX2L9XH0h9aDXYLeIXoI246DIPm10MNLsIinIMkd1GltpzcwmOrH8NsMDMtYhpXNLrC0SFdVE6Blc/XH+F/vx0mq8BK9NVNmdC/DU393Kv1uMpiwXrqFJakJD3xJ57QW/9JxR8CtrS0c1+kaRgbNSzd/XP2tp8HppTNaPt+htiVYM0Dj4DCwmbDoOVNYKzC+jrCaUlyF1XucNphxq0aR0peCu/e9C69mvVydEgVkppdwIzCxUJQcGePUMb3CaORp+OuIdhzc4u6ec62/C2JicUfAomJqJycc19kNmMODNSvAXgqzOoUpvxDmF1zMPu5Y+7cH8M10WhhfaWCZR0myV1Ui9O5pxm/ajz7UvbxfPfn+Ve7fzk6pApLTM/lw5UHmbs5ngZmIw/0bMUDhYuFOBulFPaMDD3RnziBtehbQGFXUFISlpMn4fz+f6MdsyeYA/wxNQ/D3LYL5qbN9O6fwq4gQwOpaFmbSXIX1SbHksNTvz/F7wm/M+bKMTx2zWNOO1SyLLGnspi6Yj9Ldybh7+HC+D5h3Nk9tGixkNpC2e1YT5/WE/+JRCwnErDu34zl0G4sSSexZimseQbg3JE2Rh8fTI0bYw4OLt39ExyCOSiwaksoiyolyV1UK6vdyht/vcG8A/MY3HIwr93wGi7GalwKrxqUXCykiW8DnujXhuhrmmKswBh5p2ezQtx61I4FWLYuxZqciiXPFYtbWyzGJljzXbGcPI0lKQl7evq5r9U0TAEBxRd+z34IhIQUXRQ2NmxYPbN6xUVJchfVTinFrF2z+PCfD+ka1JUP+nzg1EMly7Mh9jRvL9vHjoR02gR68kS/tlzVxIdGXi4XLWtQK9htEP9XYQXLRZBxHAwmvYJlh2HYQ/tgybBgSUzCmpRYOAT03NsqL++cXWpms97VU9bon7PdP15e1bMkYj0nyV3UmCWHl/DChhdo7tWc6f2m09izsaNDumRKKZbtSuLd5fs5nJxd9LiHi5EAL1cCvFxp5Kn/G+B53v3C2y6mWtCStdvhxD/FZRBSj4JmgOY3FBc28wo+5yVKKWxpacV9/Yklh4AWjgY6eQpstnNeZ3B3x9Q4pHAEUGHLv/C2OSQEU3AwBjfHryFQ20hyFzVqU9ImHl/9OK4mV6ZFTCO8YbijQ6oUq83On4dTSEzP5XRWAcmZ+SRn5ZOcmVd0Pz3XUuZrfRqYy03++n0XArxcaejh6hxdP0pB0o7iFv3pA4AGzboXJ3rfipV/VjYb1tOnS1/8TSr+ELCdPl3qdUY/vwt2/5gCA6u+Fn8tJ8ld1LjY1FjGrRpHWn4aU3pNoWfTno4OqVrkW22czirgdGZ+ieSfz+nCf88+djozn+wCW6nXGzTw9yhO9kU/ZXwz8HU311zXxql9xS36k7v0xxpfoyf68GHgf3m18O35+VhPniwe+1/GPAB7Zua5LzIYMAUG6om/cYh+wbfk7ZBgjP7+9ar7R5K7cIhTOad4ZNUjHEg9wIs9XuSWtrc4OiSHys63FiX985N/cmZB0YdAcmY+BTZ7qdebjVrRN4BGniU+AMr4ZuDhYqy6JHfmUHGiP7FVfyzoquJEH9Cuao5zHltWVnH3T6nyD/ptVVBwzms0V1dMwUF6l09w8DldQXopiBCMnnVnnVtJ7sJhsi3ZPPnbk2w4voEHOz7II50fqVctq8pQSpGRZy3V8k8u8YFw9sPhTHYBNnvpv1s3s6H8biFPVxqVeO6Shn2mxukVLPcu0i/MAjRqV5zog66sscJmSilsqalYTpTo8kk87/apU6XLP3h6lqr5o18QDim6EGxwrR0TvyS5C4ey2C1M/nMyPx38iaGthvLK9a9glunxVcJmV6TmFJz7TaDkN4OibwoFpGQXlLkPLzdTccIvo0vo7IdDQ0+Xc2vkZ5yAvYv1RB+3AZQd/FoWFzZrfI3DK1gqqxVrcrLe2k9MLJ4HUDQjOBFbamqp1xkbNiyn+ye4ess/XyJJ7sLhlFJ8tvMzPt76Md2Du/N+n/fxcqmGpe1EuSw2O2cKLwSfn/yLPgQKPxwy861l7sPfw6XEt4Hi6wRNXLJom/o7IceX43FiI5rdCj7N9Aux4cOhaTenrWBpz8sr7PNPLLcInP388g81Uf65AiS5C6fx86GfeWnDS7TwacGMfjMI9gi++ItEjcuz2Mq/QHz2fuFjeZZzuz18yGKA6R+GmTfRQ23HjJV0U0MONexNUpNI7M160MjHo6iryNvN5NRddUop7JmZpS7+njsPIKl6yz+XQ5K7cCp/Jv7JhDUTcDe5M73fdNr5V88FOVH9lFJk5VuLh4pmnjtUNDMjhZYp6+mas44etn9ooBVwRnmx3NaVX+zd+MN+BZrJpdR1AL1bqMQIIk83p55Ipux2bGfOlD/2P7Gw/PN5Odbg7U3Qs8/iGx1VqeNWaXLXNC0S+BAwAjOVUm+d9/xo4F3geOFD/1VKzbzQPiW51z8HUg8wbuU4sixZTO01leubXO/okEQ1s+dlkbPnV+x7FuJxdCVGazb5Ji8O+PTkb/cb2aA6ciJLcTpLv1BcVjqqzRPJlMWC5eSpUhd/vYcMwb1Ll0rts8qSu6ZpRuAA0B9IADYBo5RSe0psMxroqpR6pKIBSnKvn05mn2TcqnEcTjvMS9e9RFSbyrVeRC1kyYPDa/RJU/uXQF46uHhC24HQYRjWVhGkWM0luoEKyh09VGcmklVCRZN7Rb7vdANilVKHC3c8BxgO7Lngq4QoQ5BHEF9FfsXEtRN5aeNLJGYnMrbTWKfufxVVxOwG7QbpP9YCOPq7nuj3LYZdP2EyNSCwTT8COwzXE75bYLm7yrfaii4Ulzd0dEdCGsm1bSJZFapIcm8CxJe4nwB0L2O7WzRNuwm9lT9BKRVfxjZC4OniybR+03j1j1eZsX0GJ7JO8H/X/x9mgwyVrDdMLhDWT/8ZMhWObSwsg/Cz/mN0gdZ99VWm2g0Cd/9zXu5qMtLYtwGNfS9em/7sRLLzLxAnlygvcTg5m+SsfAqsTjSR7DJVpFtmJDBQKfVA4f27gW5KqUdLbNMQyFJK5Wua9jDwL6VU3zL29SDwIEBoaGiXuLi4qvtNRK2jlOKT7Z8wfft0rgu5jqm9p+LpUndmEopKsNsh4e/iejfp8XoFy5Y36Ym+/c3gGVAthz5/IlnpeQPFHwzlTSRrYDbSyMul7A+AEh8Qgd6uuJoqN2a+KvvcrwNeVkoNLLz/n8I34s1ytjcCKUqpC9Z+lT53cdaC2AW8svEVWvm2YnrEdII8ghwdknAGShVWsCxM9CmHiytYdhimj6f3DnFIaPbCiWRlzSAu7ibSny9rItkrw67g3utbVOrYVZncTehdLRHoo2E2AXcopXaX2CZEKZVYeDsKeEYp1eNC+5XkLkraeHwjE3+biKfZk+n9ptPWr62jQxLORCm9mNnZRJ+8T3+8WXc90YcPA99Qx8ZYjrMTyUom/2ua+xIWWLkJfVU9FHIw8AH6UMjPlVKTNU17FdislFqkadqbwDDACqQAY5VS+y60T0nu4nz7U/YzbuU4cqw5vN/nfXqEXLB9IOqz5P2FiX4hJO3UH2t8dWGiHw4NWzs2vmokk5hErZSUncTYlWM5mn6UV254hWGthzk6JOHsUg4Xt+iPb9EfC7qyuEUf0N7h9W5Kstlt2JSt0stSSnIXtVZmQSYT1k7gr8S/eKTzIzzY8UGnGYEgnFxafHEFy2N/AgoatimuYBncsUYrWCZlJ3Ew7SCxabHEpsYSmxbL4fTDvNjjRYaHDa/UfiW5i1rNYrPw8h8vs+jQIm5pcwvP93hehkqKS5OZVJzoj64vrGDZorjrpkmXKkv0Z3LP6Ak8LZaDqXoyP5R2iCxLVtE2gQ0CCfMLI8w3jMgWkVwVcFWljiXJXdR6SimmbZvG/3b8jxua3MCUXlPwMF9e0SVRT2Wf0WfF7lkIh38DuwW8mxR33TTrDoaLD03MLMjkUNohvTVe2BKPTYslJS+laBsfVx/a+LYhzDeMNn5taO3bmjDfsCpbPF6Su6gzfjrwE6/9+Rpt/NowLWIage7lz1wU4qJyU2H/Mr1FH7sKbPngEVhYqngYNL+RPGXlcPrhou6Us10rSdlJRbtpYGqgJ/HC1vjZZN7QrWG1diNKchd1yvrj63ly7ZP4uPowPWI6YX5hjg5J1AGW3FSO7ZrDwYNLiT2zh1gjxLq6EW8ycHauqtlgppVPq6Ikfjahh3iEYNBqvlCZJHdR5+w9s5dxq8aRb83ngz4f0C2km6NDErWEXdk5nnW8qCvlbEv8SPoRrHZ9kRKDZiDU1Z82FhthqScIy80iTHMjtFV/TFeM0MshmC9e7qC6SXIXddKJrBOMWzmOuMw4Xr/hdYa0GuLokIQTUUqRnJt8TldKbGosh9IPkWvNLdqusUfjUt0pLX1a4mosXEfVmg+H1+p99PuWQF4amD2g7QD9YmxYf3B1TKkMSe6izkrPT2fC2glsStrE49c8zpgrx8hQyXooPT+9aGRKyVEqGQUZRds0dGtImF9Y0QXOML8wWvu0vrQaRjYLHF1XXMEyOxlMbnrRs/CzFSyr5mJpRUhyF3Vaga2AFze8yNIjSxnZdiTPdX8Ok8E5V+wRlyfHksOhtEPF3SmFXSvJuclF23iZvUq1xFv7tsbfzf8Ce64Euw2O/VE8aSozUa9g2aq3nujbDS5VwbKqSXIXdZ5d2fl468fM3DmTnk168l6v93A3uzs6LFFJBbYCjqQfKWqJn+1aOZ51vGgbV6Nr0dDCkiNVgtyDav7bm90OxzfrXTd7FkH6MdCM0LJncWEzz6of2SXJXdQbc/fPZfJfk2nv355pEdNo1KCRo0MSF2Cz24jPjC/VEo/LiMOm9IU1TJqJFj4tilriZ7tWmng2wViB8eg1TilI3KYn+T0LIeUQoEHz64sTvU+TKjmUJHdRr/wW/xtP/f4Ufq5+zOg3g1a+rRwdUr13oen3+bZ8ADQ0mno1Pac7Jcw3jBbeLTAba+mMZKXg1J7iRJ+8V3+86bXFk6b8WlR695LcRb2z+/Ruxq0ah9Vu5cM+H9I1+KLnv6giFZp+7x54zoXNNr76CJU635V2+qCe5PcugsTt+mOD3oHuD1Vqd5LcRb2UkJnA2JVjOZ51nMk3TmZQy0GODqlOqcz0+zDfMFr7tq6y6fe1WsoRvd5N24EQ0K5Su5DkLuqt9Px0Hlv9GP+c+oeJXSYy+orRMlTyEuVZ85x2+n19V9HkLmPHRJ3j4+rDpwM+5fn1zzN1y1SOZx3nP93+45wX4hzMYrdwLONYqZZ4fGY8dqVPwD87/b5LUBenmH4vKkaSu6iTXI2uvHPTOzT2aMwXu7/gZM5J3rnpHRqYHD993BEqPP3eK5S2fm0Z3HJwUd94qFeozCGoheR/TNRZBs3AxK4TCfEM4a2/32LMr2P4uO/HNGzQ0NGhVZtLnX7fs0nPsqffi1pP+txFvbD62Gqe+f0ZGjVoxIx+M2jh08LRIV22Gpt+L5yKXFAV4jw7knfw6OpHsSkbH/f9mKsDr3Z0SBXiVNPvhcNJcheiDPEZ8YxdNZbErETeuukt+jfv7+iQilRk+r2b0Y1Wvq2cY/q9cAhJ7kKUIzUvlcdWP8b25O1M6jqJe664p0aPXyen34saI0MhhSiHn5sfnw34jOfWP8e7m98lMTuRSV0nVXnSvNTp9xGhEXVj+jFUliIAAAYBSURBVL1wCpLcRb3kZnLjvV7v8d7m9/hmzzckZifyVs+3cDO5VWp/ZU2/j02LJduSXbTN2en33YK71a/p98IhJLmLesugGXj62qdp7NGYdza9w5jl+lDJC12EvJTp90NbDZXp98JhJLmLeu+u8LsI9gjm2XXPcvfSu5nRbwaB7oEVnn7fu1lvmX4vnI5cUBWi0LZT23h09aPkWHKwKmup6fdnR6bI9HvhSHJBVYhL1DmwM7MHz+ar3V/RqEEjmX4vajU5Y4Uoobl3c1667iVHhyHEZZPvlEIIUQdJchdCiDpIkrsQQtRBktyFEKIOkuQuhBB1kCR3IYSogyS5CyFEHSTJXQgh6iCHlR/QNC0ZiKvkyxsBp6swnKoicV0aievSOWtsEteluZy4miulAi62kcOS++XQNG1zRWor1DSJ69JIXJfOWWOTuC5NTcQl3TJCCFEHSXIXQog6qLYm908dHUA5JK5LI3FdOmeNTeK6NNUeV63scxdCCHFhtbXlLoQQ4gKcNrlrmtZM07Q1mqbt1TRtt6Zpj5exjaZp2keapsVqmrZD07RrnCSu3pqmpWuatq3wp9oLhGua5qZp2t+apm0vjOuVMrZx1TTth8L36y9N01o4SVyjNU1LLvF+PVDdcZU4tlHTtK2api0u47kaf78qGJdD3i9N045qmraz8JilllFzxN9jBeOq8b/HwuP6apr2o6Zp+wrzxXXnPV+975dSyil/gBDgmsLbXsABIPy8bQYDvwAa0AP4y0ni6g0sruH3SwM8C2+bgb+AHudtMw74pPD27cAPThLXaOC/DjrPJgLflfX/5Yj3q4JxOeT9Ao4CjS7wfI3/PVYwrhr/eyw87lfAA4W3XQDfmny/nLblrpRKVEr9U3g7E9gLNDlvs+HA10r3J+CraVqIE8RV4wrfg6zCu+bCn/MvqAxHP+EAfgQitGpeybmCcTmEpmlNgSHAzHI2qfH3q4JxOasa/3t0VpqmecP/t28GL1aVYRj/PcwYxBQJFQWNMrhxYVRjMCQDLkokQWajCxeWuZEiilaiLvsDQnCRkBJRIyJTxiSjaIh/gBOCoi6GEBwsJgwVKpSpp8X5xi6HO95LeM53OLy/zT3nfO+99+Hhvs/9zvfdy0bgKIDtB7bvlMoq9aux4d5Juh0epZj1dfIScLPjfJ4ag/YRugA2pKWI05LW1aRnQNIlYAE4Z3tZv2wvAneBZxugC2BbujWdkrSqak2Jg8Be4J9lxrP41YcuyOOXgbOSZiXt6TKeqx976YL6+3EN8BvwZVpeOyJpqFRTqV+ND3dJTwHfAp/Yvlce7vKUWmaFPXT9RPEX4VeBQ8D3dWiy/bft14BhYEzSy6WSLH71oesHYMT2K8CP/DdbrgxJW4EF27OPKutyrVK/+tRVu1+JcdvrgS3Ah5I2lsZz9WMvXTn6cRBYD3xuexT4A9hXqqnUr0aHu6QVFAE6afu7LiXzQOesZRi4lVuX7XtLSxG2Z4AVkp6rWlfH+98BLgBvl4Ye+iVpEHgG+D23Ltu3bd9Pp18Ar9cgZxyYkHQDOA68KembUk0Ov3rqyuQXtm+lxwXgJDBWKsnSj710ZerHeWC+4y51iiLsyzWV+dXYcE9rm0eBa7Y/W6ZsGng37Tq/Ady1/UtuXZJeXFqblTRG4fPtinU9L2llOn4S2ARcL5VNA7vS8XbgvNPOTk5dpXXGCYp9jEqxvd/2sO0Ris3S87Z3lspq96sfXTn8kjQk6emlY2AzcKVUlqMfe+rK0Y+2fwVuSlqbLr0FXC2VVerX4ON6oQoYB94BLqf1WoADwGoA24eBGYod5zngT2B3Q3RtBz6QtAj8BeyoOhQofsXzlaQBig/vCdunJH0KXLQ9TfGl9LWkOYoZ6I6KNfWr62NJE8Bi0vVeDbq60gC/+tGVw68XgJMpIweBY7bPSHofsvZjP7py9CPAR8CkpCeAn4HddfoV/1ANgiBoIY1dlgmCIAj+PxHuQRAELSTCPQiCoIVEuAdBELSQCPcgCIIWEuEeBEHQQiLcgyAIWkiEexAEQQv5F72fAXn94sJXAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "n_comp = range(70, 101, 10)\n",
    "K_pca = range(2, 8, 2)\n",
    "\n",
    "scs_pca=[]\n",
    "chis_pca=[]\n",
    "\n",
    "for n in n_comp:\n",
    "    pca = PCA(n_components=n)\n",
    "    #print(\"PCA begin with n_components: {}\".format(n));\n",
    "    pca.fit(train)\n",
    "    \n",
    "    train_pca = pca.transform(train)\n",
    "    \n",
    "    for K_0 in K_pca:\n",
    "        chi, sc = K_cluster_analysis(K_0, train_pca)\n",
    "        chis_pca.append(chi)\n",
    "        scs_pca.append(sc)\n",
    "    \n",
    "    plt.plot(K_pca,np.array(scs_pca),label=\"n_components=%d\"%(n,))\n",
    "    scs_pca.clear()\n",
    "    \n",
    "leg = plt.legend(loc='best', ncol=2, mode=\"expand\", shadow=True, fancybox=True)\n",
    "leg.get_frame().set_alpha(0.5)\n",
    "plt.show() \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "由图可见，降维到80，并且K选择为2时，得分最高。图片显示的是轮廓系数。根据打印的CH_score，可见最优点是几乎一致的。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "这次作业k-means方法运行了几次，但每次结果都有所不同。这可能是符合k-means的工作原理——即受初始随机点的影响最终聚类的结果。但k-means的聚类是适合球型分布的，我很想把散点图画出来看一看，但是不会画。因为这次没有标签集，老师的示例代码又有奇怪的错误，时间不足，所以以后有机会再尝试。如果样本不是高斯分布的，最好再尝试使用DBSCAN聚类，看一下效果。"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
